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IV 


ABSTRACT 


In  this  thesis,  the  impacts  of  transforming  the  coordinate  system  of  an  existing 
x-z  mesoscale  model  to  x-oy  are  analyzed  and  discussed  as  they  were  observed  in 
three  test  cases.  The  three  test  cases  analyzed  are:  A  rising  thermal  bubble,  a 
linear  hydrostatic  mountain,  and  a  linear  nonhydrostatic  mountain.  The  methods 
are  outlined  for  the  transformation  of  two  sets  (set  1,  the  non-conservative  form 
using  Exner  pressure,  momentum,  and  potential  temperature;  and  set  2,  the  non¬ 
conservative  form  using  density,  momentum,  and  potential  temperature)  of  the  x-z 
Navier-Stokes  equations  to  x-oy  and  their  spatial  (Continuous  Galerkin)  and  temporal 
(Runge-Kutta  35)  discretization  methods  are  shown  in  detail.  For  all  three  test  cases 
evaluated,  the  x-oy  models  performed  worse  than  their  x-z  counterparts,  yielding 
higher  RMS  errors,  which  were  observed  predominantly  in  intensity  values  and  not 
in  placement  of  steady  state  features.  Since  the  models  did  converge  to  a  fairly 
representative  steady-state  solution  the  results  found  by  this  project  are  promising, 
even  though  they  did  indicate  that  x-oy  coordinates  are  not  as  accurate  or  efficient 
as  x-z  coordinates.  With  further  fine-tuning  of  the  model  environment,  these  issues 
could  be  made  minimal  enough  to  warrant  their  utility  with  semi-implicit  methods. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

II.  BACKGROUND .  5 

A.  GOVERNING  EQUATIONS .  5 

1.  Equation  Set  1:  Navier-Stokes  Equations  with  Exner 

Pressure .  5 

2.  Equation  Set  2:  Navier-Stokes  Equations  with  Density  .  6 

B.  X-Z  TO  X-az  COORDINATE  SYSTEM  TRANSFORM  ....  7 

1.  Gal- Chen  and  Somerville .  7 

2.  Basic  Transformation  Machinery .  9 

3.  Transformation  Functions  .  12 

C.  SPATIAL  DISCRETIZATION:  CONTINUOUS  GALERKIN  .  .  14 

1.  Equation  Set  1 .  18 

2.  Equation  Set  2 .  18 

D.  TEMPORAL  DISCRETIZATION  RUNGE-KUTTA  35 .  19 

III.  APPLIED  COORDINATE  TRANSFORMS .  21 

A.  EQUATION  SET  1 .  21 

1.  Perturbation  Method .  21 

2.  Transform .  23 

3.  Decomposition .  26 

4.  Application  of  the  Galerkin  Statement  .  27 

B.  EQUATION  SET  2 .  29 

1.  Perturbation  Method .  29 

2.  Transform .  31 

3.  Decomposition .  34 

4.  Application  of  the  Galerkin  Statement  .  35 

vii 


IV.  TEST  CASES .  37 

A.  OVERVIEW .  37 

B.  CASE  1:  RISING  THERMAL  BUBBLE .  37 

C.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN .  41 

D.  CASE  3:  LINEAR  NON-HYDROSTATIC  MOUNTAIN .  45 

V.  RESULTS .  47 

A.  CASE  1:  RISING  THERMAL  BUBBLE .  47 

1.  Accuracy  and  Comparison .  47 

2.  Conclusions .  49 

B.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN .  53 

1.  Accuracy  and  Comparison .  53 

2.  Conclusions .  56 

C.  CASE  3:  LINEAR  NON-HYDROSTATIC  MOUNTAIN .  59 

1.  Accuracy  and  Comparison .  59 

2.  Conclusions .  63 

VI.  CONCLUSIONS  AND  RECOMMENDATIONS .  67 

APPENDIX.  COEFFICIENTS  FOR  RK35  METHOD .  69 

LIST  OF  REFERENCES  .  71 

INITIAL  DISTRIBUTION  LIST  .  75 


viii 


LIST  OF  FIGURES 


1.  x-z  to  x-cr2  coordinate  transformation:  (i)  traditional  x-z  coordinates, 

(ii)  x-cr2  coordinates,  and  (iii)  x-cr2  coordinates  mapped  back  to  x-z  space.  9 

2.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  per¬ 

turbations  using  unmodified  CGI  source  codes  [1]  after  700  s  for  reso¬ 
lutions:  (a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.05  to  0.525  K  with  an  interval  of 
0.025  K .  39 

3.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  per¬ 
turbations  using  unmodified  CG2  source  codes  [1]  after  700  s  for  reso¬ 
lutions:  (a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.05  to  0.525  K  with  an  interval  of 


0.025  K .  40 

4.  Case  2:  Linear  Hydrostatic  Mountain.  The  single  mountain  peak  zsurf, 

(a),  and  the  associated  slope  of  the  peak  dZg^rf ,  (b) .  43 

5.  Case  2:  Linear  Hydrostatic  Mountain.  Resulting  horizontal  velocity  (a) 


and  vertical  velocity  (b)  using  unmodified  (i)  CGI  and  (ii)  CG2  source 
codes  [1]  after  10  h  for  the  resolution  of  1200  nr  (in  x)  and  240  m  (in 
z)  (colored  lines)  plotted  with  the  analytic  solution  (black  lines).  All 
cases  were  run  using  10th  order  polynomials,  with  contours  from  -0.025 
to  0.025  ms-1  with  an  interval  of  0.005  ms-1,  (a),  and  -0.005  to  0.005 
ms-1  with  an  interval  of  0.0005  ms-1,  (b) .  44 


IX 


6.  Case  3:  Linear  Non- Hydrostatic  Mountain.  Resulting  horizontal  veloc¬ 
ity  (a)  and  vertical  velocity  (b)  using  unmodified  (i)  CGI  and  (ii)  CG2 
source  codes  [1]  after  5  h  for  the  resolution  of  360  m  (in  x)  and  300 
m  (in  z)  (colored  lines)  plotted  with  the  analytic  solution  (black  lines). 

All  cases  were  run  using  10th  order  polynomials,  with  contours  from 
-0.025  to  0.025  ms-1  with  an  interval  of  0.005  ms'1,  (a),  and  -0.005  to 
0.005  ms'1  with  an  interval  of  0.0005  ms-1,  (b) .  46 

7.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  per¬ 

turbations  using  modified  CGI  source  codes  [1]  after  700  s  for  resolu¬ 
tions:  (a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.05  to  0.525  K  with  an  interval  of 
0.025  K .  50 

8.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  per¬ 

turbations  using  modified  CG2  source  codes  [1]  after  700  s  for  resolu¬ 
tions:  (a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.05  to  0.525  K  with  an  interval  of 
0.025  K .  51 

9.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  per¬ 

turbations  for  all  four  models  along  the  vertical  axis  (x  =  500  m)  after 
700  s  for  resolutions:  (a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run 
using  10th  order  polynomials .  52 

10.  Case  2:  Linear  Hydrostatic  Mountain.  Resulting  horizontal  velocity  (a) 
and  vertical  velocity  (b)  using  modified  (i)  CGI  and  (ii)  CG2  source 
codes  [1]  after  10  h  for  the  resolution  of  1200  m  (in  x)  and  240  m  (in 
z)  (colored  lines)  plotted  with  the  analytic  solution  (black  lines).  All 
cases  were  run  using  10th  order  polynomials,  with  contours  from  -0.025 
to  0.025  ms'1  with  an  interval  of  0.005  ms-1,  (a),  and  -0.005  to  0.005 

ms'1  with  an  interval  of  0.0005  ms”1,  (b) .  57 


x 


11.  Case  2:  Linear  Hydrostatic  Mountain.  Normalized  momentum  flux  for 
the  resolution  of  1200  m  (in  x)  and  240  m  (in  z),  at  times  2  h,  4  h,  6 
h,  8  h,  and  10  h  for  the  four  model  runs:  (i)  CGI  x-z,  (ii)  CGI  x-crz, 

(iii)  CG2  x-z,  and  (iv)  CG2  x-az .  58 

12.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Resulting  horizontal  veloc¬ 
ity  (a)  and  vertical  velocity  (b)  using  both  modified,  (ii),  and  unmodi¬ 
fied,  (i),  CGI  source  codes  [1]  after  5  h  for  the  resolution  of  360  m  (in 
x)  and  300  m  (in  z)  (colored  lines)  plotted  with  the  analytic  solution 
(black  lines).  All  cases  were  run  using  10th  order  polynomials,  with 
contours  from  -0.025  to  0.025  ms-1  with  an  interval  of  0.005  ms-1,  (a), 

and  -0.005  to  0.005  ms-1  with  an  interval  of  0.0005  ms-1,  (b) .  60 

13.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Resulting  horizontal  veloc¬ 
ity  (a)  and  vertical  velocity  (b)  using  both  modified,  (ii),  and  unmodi¬ 
fied,  (i),  CG2  source  codes  [1]  after  5  h  for  the  resolution  of  360  m  (in 
x)  and  300  m  (in  z)  (colored  lines)  plotted  with  the  analytic  solution 
(black  lines).  All  cases  were  run  using  10th  order  polynomials,  with 
contours  from  -0.025  to  0.025  ms-1  with  an  interval  of  0.005  ms-1,  (a), 

and  -0.005  to  0.005  ms-1  with  an  interval  of  0.0005  ms-1,  (b) .  61 

14.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Normalized  momentum 
flux  for  the  resolution  of  360  m  (in  x)  and  300  m  (in  z),  at  times  1  h, 

2  h,  3  h,  4  h,  and  5  h  for  the  four  model  runs:  (i)  CGI  x-z,  (ii)  CGI 
x-crz,  (iii)  CG2  x-z,  and  (iv)  CG2  x-crz .  65 


xi 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


LIST  OF  TABLES 


I.  Case  1:  Rising  Thermal  Bubble.  Comparison  of  modified  and  unmod¬ 
ified  CGI  and  CG2  using  5  m  resolution  after  700s .  48 

II.  Case  1:  Rising  Thermal  Bubble  RMSE.  Root-mean-squared  errors  for 

the  four  primary  state  variables,  for  the  modified  codes  using  the  un¬ 
modified  code  solutions  as  the  analytic  solutions,  after  700  s  using  5  m 
resolution  and  10th  order  polynomials .  49 

III.  Case  2:  Linear  Hydrostatic  Mountain.  Comparison  of  modified  and  un¬ 
modified  CGI  and  CG2  using  1200  m  (in  x)  and  240  m  (in  z)  resolution 

after  10  hours .  54 

IV.  Case  2:  Linear  Hydrostatic  Mountain  RMSE.  Root- mean-squared  er¬ 

rors  for  the  four  primary  state  variables,  for  both  the  modified  and 
unmodified  codes,  after  10  h  using  1200  m  (in  x)  and  240  m  (in  z) 
resolution  and  10th  order  polynomials .  55 

V.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Comparison  of  modified 

and  unmodified  CGI  and  CG2  using  360  m  (in  x)  and  300  m  (in  z) 
resolution  after  5  h .  62 

VI.  Case  3:  Linear  Non-Hydrostatic  Mountain  RMSE.  Root-mean-squared 

errors  for  the  four  primary  state  variables,  for  both  the  modified  and 
unmodified  codes,  after  5  h  using  360  m  (  in  x)  and  300  m  (in  z) 
resolution  and  10th  order  polynomials .  64 


xiii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


ACKNOWLEDGMENTS 


I  would  like  to  thank  my  advisors,  Francis  X.  Giraldo  and  Major  Tony  Eckel, 
without  whom  this  thesis  would  not  have  been  possible.  I  would  also  like  to  thank 
the  Naval  Postgraduate  School,  the  Meteorology  and  Mathematics  departments  at 
NPS,  The  Air  Force  Institute  of  Technology,  The  United  States  Air  Force,  and  The 
United  States  Navy. 


xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xvi 


I. 


INTRODUCTION 


Numerical  Weather  Prediction  (NWP)  models  are  the  work  horse  of  mod¬ 
ern  atmospheric  condition  forecasting  and  as  such  have  been  the  focus  of  numerous 
studies,  aiming  to  develop  more  accurate  models  with  longer  deterministic  forecast 
periods.  There  are  countless  areas  in  which  to  focus  this  research  (spatial  and  tem¬ 
poral  discretization  methods,  non-reflecting  boundary  conditions,  data  assimilation, 
prognostic  equations,  physical  parameterizations,  etc....),  and  for  this  project  the  fo¬ 
cus  will  be  on  the  transformation  of  the  prognostic  or  governing  equations  from  x-z 
coordinates  to  terrain  following  x-az  coordinates  using  a  specific  spatial  (continu¬ 
ous  Galerkin)  and  temporal  (Runge-Kutta  35)  discretization  method.  This  thesis  is 
necessary  in  order  for  future  research  to  be  able  to  evaluate  and  compare  various 
coordinates  systems  while  varying  the  temporal  discretization  methods,  allowing  for 
larger  time  steps  while  maintaining  stability  (i.e.,  semi-implicit  methods). 

Various  mature  mesoscale  models  such  the  Coupled  Ocean/ Atmosphere  Meso- 
scale  Prediction  System  (COAMPS)  [2]  and  the  Weather  Research  and  Forecasting 
(WRF)  modeling  system  [3]  use  a  variation  of  the  x-cr  coordinate  transformation.  By 
studying  the  works  of  Gal-Chen  and  Somerville  [4]  and  analyzing  the  transforma¬ 
tion  methods  used  in  COAMPS  and  WRF,  a  similar  x-crz  coordinate  transformation 
was  applied  to  the  governing  equations  of  interest  employing  continuous  Galerkin 
techniques. 

The  governing  equations  selected  for  this  thesis  consist  of  two  forms  of  the 
Navier-Stokes  equations.  The  Navier-Stokes  equations,  along  with  their  variations, 
form  the  most  widely  used  and  accepted  sets  of  equations  for  numerically  resolving 
atmospheric  flow.  The  first  specific  formulations  of  the  equation  sets  selected  (set 
1)  was  the  non-conservative  form  using  Exner  pressure,  momentum,  and  potential 
temperature,  which  is  used  in  the  operational  NWP  model  COAMPS  [2],  The  second 
formulation  chosen  (set  2)  was  the  non- conservative  form  using  density,  momentum, 
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and  potential  temperature.  The  operational  NWP  model  WRF  [3]  also  uses  set  2, 
but  in  a  conservative  form.  Building  on  the  work  of  Giraldo  [1]  who  by  implementing 
continuous  Galerkin  techniques,  developed  a  2-D  (x-z  slice)  mesoscale  model  using 
Non-Hydrostatic  Equations  (Euler  and  Navier-Stokes  Equations),  the  original  con¬ 
struct  was  transformed  from  x-z  coordinates  to  x-az  coordinates  to  test  the  impacts 
on  resolving  atmospheric  motion  in  a  continuous  Galerkin  (CG)  framework.  Cur¬ 
rently,  most  operational  non-hydrostatic  models  use  finite  difference  (FD)  methods 
(i.e.,  structured  grids),  which  then  rely  on  x-a  in  order  to  resolve  atmospheric  flow 
in  the  presence  of  terrain.  CG  methods,  on  the  other  hand,  can  use  various  types 
of  grids  (unstructured  grids,  x-z  grids,  x-az  grids,  etc....).  For  this  reason,  the  CG 
method  is  well  positioned  to  judge  the  effects  of  the  coordinate  system  on  the  solution 
accuracy  and  efficiency,  which  is  the  goal  of  this  thesis. 

After  the  modifications  were  made  to  the  model,  three  test  cases  were  run:  ris¬ 
ing  thermal  bubble,  linear  hydrostatic  mountain,  and  linear  non-hydrostatic  moun¬ 
tain.  The  numerical  solutions  were  either  evaluated  against  other  model  solutions 
(case  1)  or  the  analytic  approximations  (case  2  and  case  3)  using  root  mean  squared 
error  and  normalized  momentum  flux.  The  resultant  data  was  also  compared  to 
the  unmodified  solutions.  Additionally,  the  initial  conditions  for  the  test  cases  were 
pre-defined  for  each  case,  maintaining  uniform  initial  conditions  from  which  both 
coordinate  systems  numerical  solutions  can  be  compared. 

With  x-az  coordinates  that  are  proven  to  function  properly  using  fully  explicit 
time  integration,  future  research  will  be  able  to  evaluate  x-az  coordinates  using  semi- 
implicit  methods.  NWP  models  are  already  taking  advantage  of  semi-implicit  time 
integration  methods,  which  optimize  the  horizontal  and  vertical  resolutions  and  their 
associated  time  sets  while  maintaining  stability.  The  2-D  semi-implicit  method  (x- 
z  formulation)  can  use  very  large  time-steps  sizes  since  the  Courant-Friedrichs-Lewy 
(CFL)  condition  is  no  longer  constrained  by  acoustic  and  gravity  waves;  the  penalty  is 
that  a  global  2-D  implicit  problem  must  be  solved.  In  contrast,  the  1-D  semi-implicit 


2 


x-az  formulation  can  use  large  time-steps,  but  it  must  adhere  to  the  CFL  condition 
due  to  gravity  and  acoustic  waves  in  the  horizontal;  the  advantage  is  that  only  a  1-D 
matrix  problem  needs  to  be  solved.  A  model  using  x-az  can  exploit  semi-implicit 
methods  along  the  vertical  (a)  direction,  which  is  the  inspiration  of  this  thesis  topic. 
Both  COAMPS  and  WRF  are  currently  using  semi-implicit  methods  (only  along  a), 
but  are  constrained  by  FD  spatial  discretization  methods. 

Though  only  explicit  time  integration  was  used  in  this  thesis,  it  is  now  possible 
to  observe  if  the  implementation  of  x-az  coordinates  significantly  improves  or  dimin¬ 
ishes  the  solution  of  the  Navier-Stokes  equations  over  x-z  coordinates  when  using  a 
CG  framework.  Giraldo  [5]  has  already  developed  semi-implicit  methods  for  the  2-D 
model  in  the  x-z  framework  and  tested  their  impacts.  With  the  model  coordinates 
transformed  to  x-crz,  future  research  will  be  able  to  extent  the  time  integration  to 
semi-implicit  methods  for  x-cr2  and  compare  the  results  with  the  semi-implicit  results 
derived  using  x-z  coordinates.  The  relevant  governing  equations  for  this  project,  the 
coordinate  transformation  theory,  and  the  discretization  methods  are  discussed  in 
Chapter  II.  The  application  of  the  coordinate  transforms  are  discussed  in  Chapter 
III.  The  three  test  cases  are  explained  and  outlined  in  Chapter  IV.  A  discussion  and 
interpretation  of  the  resulting  impacts  from  the  transformation  is  in  Chapter  V.  The 
conclusions  found  and  recommendations  are  presented  in  Chapter  VI. 
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II. 


BACKGROUND 


A.  GOVERNING  EQUATIONS 


1.  Equation  Set  1:  Navier-Stokes  Equations  with 
Exner  Pressure 

The  first  set  of  equations  chosen  was  the  Navier-Stokes  equations  that  uses 
Exner  pressure,  of  which  there  has  been  extensive  amounts  of  documented  work  (i.e. 
this  is  the  formulation  used  in  COAMPS)  for  comparison.  This  set  can  only  be  writ¬ 
ten  in  non-conservative  form  and  consists  of  a  system  of  three  equations.  The  first 
equation  is  the  pressure  tendency  equation: 


8tt  _  _  R  _ 

— — b  u  ■  V 7r  H - 7rV  ■  u  —  0 

at  c„ 


(2.1) 


71  = 


JL 

Cp 


where  tt  is  Exner  pressure,  u  represents  the  velocity  field  (u,  w),  R  is  the  gas  con¬ 
stant,  cv  is  the  specihc  heat  for  constant  volume,  cp  is  the  specihc  heat  for  constant 
pressure,  P  is  pressure,  and  P0  is  pressure  at  the  surface.  The  second  equation  is  the 
momentum  equation: 


du 

dt 


+  u  ■  X7u  +  CpdV Ti 


— gk  +  /j\72u 


(2.2) 


where  6  is  potential  temperature,  g  is  the  gravitational  constant,  k  is  a  vector  (0, 1)T, 
and  yU  is  the  dynamic  viscosity.  The  third  equation  is  the  thermodynamic  energy 
equation: 


5 


—  +  u-  ve  =  pV2e 

at 


(2.3) 


For  the  scope  of  this  project,  only  inviscid  flow  will  be  considered  (i.e.  n  —  0)  and  the 
equations  further  reduce  to  the  Euler  equations  which  will  be  used  in  the  following 
sections. 


2.  Equation  Set  2:  Navier-Stokes  Equations  with  Den¬ 
sity 

The  second  set  of  equations  chosen  was  a  version  of  the  Navier-Stokes  equa¬ 
tions  that  is  now  used  in  contemporary  NWP  models  (i.e.  this  is  the  formulation 
used  by  WRF)  and  uses  density,  momentum,  and  potential  temperature  as  the  pri¬ 
mary  state  variables.  Unlike  set  1,  this  set  can  be  written  in  both  conservative  and 
non- conservative  form.  Though  neither  form  can  conserve  energy,  using  the  non¬ 
conservative  form  can  still  conserve  mass  and  more  sophisticated  time  integration 
strategies  can  be  used  [6].  Thus  for  this  thesis  the  non- conservative  form  will  be  used, 
which  consists  of  a  system  of  three  equations.  The  Erst  equation  is  the  mass  equation: 

d  +  V  •  (pit)  =  0  (2.4) 


where  p  is  density.  The  second  equation  is  the  momentum  equation: 


du  1  — * 

—  +  u  ■  V«  +  -  VP  =  —gk  +  uV’m 
at  p 


P 


Po 


m 


(2.5) 
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where  P  is  pressure,  Pq  is  pressure  at  the  surface,  and  7  is  — .  The  third  equation  is 

Cy 

the  thermodynamic  energy  equation: 

r)0 

-wr  +  u-VO  =  pV29  (2.6) 


Similar  to  set  1,  only  inviscid  flow  will  be  considered  (i.e.  p  =  0)  and  the  equations 
further  reduce  to  the  Euler  equations  which  will  be  used  in  the  following  sections. 

B.  X-Z  TO  X-az  COORDINATE  SYSTEM  TRANSFORM 

1.  Gal-Chen  and  Somerville 

In  1975,  Gal-Chen  and  Somerville  took  the  anelastic  approximation  of  the 
Navier-Stokes  Equation  (in  the  cartesian  form)  and  transformed  the  coordinated  sys¬ 
tem  to  sigma-z  coordinates  [4],  An  expanded  set  of  prognostic  equations  was  used 
for  Gal-Chen  and  Somerville’s  derivation  and  only  the  first  three  equations  and  the 
resulting  transform  will  be  used  for  comparison  in  this  project.  The  first  equation 
was  the  continuity  equation: 


G PouJ),j=  0. 


(2.7) 


where  p0  is  density,  u  is  the  velocity  components  (u,  v,  w  )T ,  and  j  is  the  derivative 
operator  where: 


d 

(Pou3),j=  g-{PoU3) 
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The  second  equation  was  the  momentum  equation: 


(po«*)  +  (poll'll-’), i  =  -(6 vp'),j  +5,3p'g  +  Tv,j . 


(2.8) 


where  SlJ  is  the  Kronecker  delta,  p'  is  the  density  perturbation,  p'  is  the  pressure 
perturbation,  and  r*-7  is  the  eddy  viscosity.  The  third  equation  used  was  the  thermo¬ 
dynamic  energy  equation: 


St 


(po0')  +  (po0'uj),j=H*  [ 


(2.9) 


where  6'  is  the  perturbation  in  potential  temperature  and  H i  is  the  eddy  diffusion. 

Using  equations  2.7  -  2.9,  Gal-Chen  and  Somerville  derived  the  following  set 
of  transformations  [4]: 


x  =  x,  y  =  y, 


_  _  H(z  -  zs) 
Z~  ( H-z8 ) 


dz_  _  dzs  z-H  dz_  _  dzs  z  —  H  dz_  _  H 

dx  dx  H  —  Zg  dy  dy  H  —  zs'  dz  H  —  zs 


u  ’ 

V 

= 

,  w  , 

1, 

0, 


0, 

1, 


dzs  z—H  dzs  z—H 


o 

0 

z-H 


V 

\W  J 


where  the  variables  with  the _  represent  the  variables  that  have  been  transformed,  H 
is  the  height  at  the  top  of  model  space,  and  zs  is  the  height  at  the  surface  of  the 
model. 


i)  ii)  iii) 

Figure  1.  x-z  to  x-crz  coordinate  transformation:  (i)  traditional  x-z  coordinates,  (ii) 
x-az  coordinates,  and  (iii)  x-az  coordinates  mapped  back  to  x-z  space. 


These  transformation  functions  are  used  to  convert  traditional  x-z  (see  figure 
l.i)  to  x-az  coordinates  (see  figure  l.ii),  which  when  mapped  back  to  x-z  space  are 
terrain  following  (see  figure  l.iii).  The  derived  inverse  transformations  are: 


z(H-zs) 

x  =  x,  y  =  y,  z  =  [ - — - J  +  z. 


u 

V 

= 

,  w  . 

0, 

dzs  z—H 

’  dx  H  ’ 


o, 

i. 


0 

0 


dz3  z—H  H—zs 
'  dy  H  ’  H  J 


V 

\W  J 


2.  Basic  Transformation  Machinery 

This  section  will  outline  the  basic  equations  used  to  transform  the  x-z  coordi¬ 
nates  of  the  Navier-Stokes  Equations  to  x-az,  which  are  similar  to  the  transformations 
derived  by  Gal-Chen  and  Somerville  in  the  previous  section,  but  in  two  dimensions. 
The  first  concept  used  was  the  total  differential: 


dx  dx 

dx  =  —dx  +  -7—dz 
dx  dz 
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da 


daz 

dx 


dx  + 


daz 

dz 


dz 


where  the' notation  indicates  the  transformed  variable.  The  two  total  derivatives  can 
then  be  written  as  a  system  of  equations: 


^  dx 
y  da  z 


/ 

dx 

dx 

\ 

dx 

dz 

da  z 

da  z 

V 

dx 

dz 

) 

dx 
y  dz 


(2.10) 


Using  vector  notation  (^)  the  above  system  can  be  simplified  to: 


dx  =  Jdx 


(2.11) 


where  J  is  the  Jacobian  of  the  transformation.  Considering  that  the  velocity  u  is 
defined  by  the  transformation  can  then  be  manipulated  to  become: 


u  =  Ju 


(2.12) 


Using  basic  linear  algebra  and  assuming  that  J  is  invertible,  which  is  true  for  an  affine 
(one-to-one)  mapping,  the  inverse  transform  for  velocity  can  be  written  as: 


u  =  J  1u 


(2.13) 


The  next  major  mathematical  concept  derived  for  the  transformation  was  the  gradi¬ 
ent  operator  (V)  which,  as  a  system  of  equations,  can  be  written  as: 
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d  d  dx  d  daz 

dx  dx  dx  daz  dx 

d  d  dx  d  daz 

dz  dx  dz  daz  dz 


which  when  put  into  matrix  form  becomes: 


(a) 

dx  daz  ^ 

(  d_  ) 

dx 

dx  dx 

dx 

d 

dx  daz 

d 

\  dz  ) 

\  dz  dz  J 

\  daz  J 

Similar  to  the  total  derivative,  the  gradient,  using  vector  notation,  was  the  defined  as: 


V  =  JTV 


(2.14) 


The  last  major  concept  used  for  the  coordinate  transformation  is  the  linear  algebra 
identity: 


(. AB)t  =  BtAt 


(2.15) 


where  A  e  7 Zmxn  and  B  e  lZnxl.  The  matrices  were  multiplied  together  and  then 
transposed  and  also  rearranged,  transposed,  and  then  multiplied  to  show  equality: 


(ABf  = 


l  T 


Ol,l  ••• 


\  4  •••  am’n  J  y  ^n,l  •••  J 


bi,i  ...  &i,, 
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Q'l, 1^1,1  +  •••  +  Ol,n&n,l  •••  al,\bl,l  +  •••  +  Ol  ,n.bn,l 

am,  1&1,1  +  •••  +  am,nbn,  1  •••  am,lbl,l  +  •••  +  am,nbn,l 


Ol,  1^1,1  + 

...  +  CLl,nbn,i 

Om,  1^1,1  +  ■ 

•  •  Clm,nbn,l 

1 - 

+ 

. . .  H-  cti  nbn  i 

0>m,lbl,l  +  ■ 

. .  +  CLm,nbn,l 

(  b 
°1,1 

b  \ 
0n,  1 

( 

...  CLm\ 

BtAt 

y  bij 

bn,i  j 

y  0  1  ,n 

...  0"m,n  J 

Ol, 1^1,1  +  •• 

•  +  Ol,n^n,l 

0>m,  it 

'1,1  +  - 

•  H- 

ai,i&i,z  +  •• 

•  +  Ol  ,nbn,l 

it 

\l  +  - 

•  “1“  QJm,nbn,l 

3.  Transformation  Functions 

For  the  transformations  from  x-z  coordinates  to  x-az,  the  x  and  z  coordinates 
from  Gal-Chen  and  Somerville  [4]  were  rewritten  using  new  notation  which  will  be 
seen  through  the  remainder  of  the  thesis: 


x  =  x,  <7  z 


H(z  -  zs) 
(H  -  zs) 


(2.16) 


where  x  is  now  x  and  z  is  a z.  Using  Eqs.  (2.16)  their  derivatives  were  evaluated, 
yielding: 


daz  dzs  crz  —  H  daz  H 

dx  dx  H  —  zs’  dz  H  —  zs 


(2.17) 
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which  matched  the  derivatives  found  by  Gal-Chen  and  Somerville  [4],  The  derivatives 
were  then  applied  to  various  forms  of  the  Jacobian  from  Eq.  (2.10),  which  yielded: 


J 


1,  0 

dzs  az-H  H 

\  dx  H-zs  ’  H — Zs 


(2.18) 


JT 


dzs  az—H 
dx  H—zs 


\ 


o, 


H 

H-Zs 


(2.19) 


J"1 


/ 


1, 


0 


l  _ dza  u  z  H  H—zs 

\  dx  H  ’  H 


(2.20) 


(J~T 


dzs  crz~H 
dx  H 


H—zs 

H 


(2.21) 


Applying  Eq.  (2.18)  to  Eq.  (2.12)  gives  the  transformation  functions  for  the  velocity 
held: 


a  ) 

/ 

1, 

0 

w  J 

dzs 

<Tz-H 

H 

\  dx 

H—zs  ’ 

H-Zs 

y  w 


(2.22) 


Using  Eq.  (2.20)  and  some  algebraic  manipulation  the  inverse  transformation  func¬ 
tions  can  be  constructed: 


<vz(H- 

X  =  x,  y  =  y,  z=  [ - — 


M 

( 

i, 

0 

“ ) 

dzs  <Tz—H 

H—zs 

l ' 

dx  H  > 

H 

/ 

V 


(2.23) 


(2.24) 
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C.  SPATIAL  DISCRETIZATION:  CONTINUOUS  GALERKIN 

Continuous  Galerkin  (CG)  methods  are  the  general  family  of  methods  that 
include:  finite  element  (FE),  spectral  element  (SE),  and  spectral  methods.  CG  meth¬ 
ods  take  complex  geometries  (elements)  from  physical  space  to  computational  space, 
using  continuous  basis  functions  that  are  used  to  approximate  the  solution  to  a  given 
partial  differential  equation  (PDE).  To  construct  the  problem  relevant  to  this  project, 
the  generalized  2-D  hyperbolic-elliptic  PDE  was  first  considered  [7]: 

dq  2 

—  +  u  ■  Vq  —  vv  q 


where  q  =  q(x,t),  u  =  u(x),  x  =  ( x,z)T ,  and  v  is  the  viscosity  coefficient.  Using 
Galerkin  machinery,  q  and  u  were  approximated  using  basis  function  expansions: 


mn 

qN(x,t)  =  '52'il>j(x)qj(t)  (2.25) 

3= 1 

Mjv 

uN(x,t)  =  ^^j(x)uj(t)  (2.26) 

3= 1 


where  Mjv  represents  the  number  of  points  inside  the  quadrilateral  elements  (Mn  = 
(N  +  l)2),  N  is  the  order  of  the  polynomial  approximation,  and  tpj  is  the  Lagrange 
polynomial  basis  functions.  The  approximations  for  qN  and  Un  were  then  substituted 
into  the  PDE,  multiplied  by  a  test  function,  i/jj,  and  integrated  across  the  global 
domain,  Q,  to  get  (weak  integral  form): 


/  i;I^dn+  [  ^(u-  VqN)dn  =  V  [  tj)I'V2qNdn  VT  G  H1  (2.27) 
it!  at  Jn  Jn 
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Instead  of  solving  the  global  problem  directly,  2-D  local  basis  functions  were  con¬ 
structed.  The  2-D  local  basis  functions  are  defined  as: 


=  MO  ®  MO 

ji  k  =  0,1, N  i  —  1,2, ...,  (N  +  l)2 

where  h  is  a  1-D  local  basis  function  and  <S>  is  the  tensor/outer  product  of  the  1- 
D  local  basis  functions.  The  1-D  Lagrange  polynomial  local  basis  function,  using 
Legendre-Gauss-Lobatto  interpolation  points,  was  defined  by: 

*'(f)=  n  (£ |) 

I  =  o 

l*j 

Additionally,  in  order  to  construct  the  basis  functions  and  transition  between  physi¬ 
cal  space  (x,  z )  and  computational  space  (£,  rf)  requires  knowledge  of  the  metric  terms: 

d£  1  dz  d£  —1  dx 
dx  \J\dr)  ’  dz  \J \  drj 

dr)  —  1  dz  dr)  1  dx 

dx  \J\  <9£  ’  dz  \J\  d£ 

dx  dz  dx  dz 
d£  dr)  dr)  d£ 

With  Eq.  (2.27)  in  local  form  it  can  now  be  converted  using  integration  by  parts 
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(IBP)  giving: 


[  ^V2qNdne  =  I  V  •  (^iVqN)dQe  -  I  V ipi  •  V (qN)dtte  (2.28) 


then  using  the  divergence  theorem  on  Eq.  (2.28)  yields: 


[  i/),V2qNdIle  =  [  n  ■  qN)dV e  -  [  V^i  •  V(qN)dQe  (2.29) 

Jre 


The  result  from  Eq.  (2.29)  is  then  substituted  back  into  the  original  PDE  seen  in  Eq. 
(2.27),  which  produces: 


VqN)dtte 


(^iVqN)dYe  -  v  f  •  V (qN)dn 
J  Qe 


e 


Substituting  the  summation  approximation  for  qN  and  (Eqs.  2.25  and  2.26)  yields: 


(mn  \ 

I  Ev^  I  dne 


(mn  \ 

( £  v '  yp  j 


The  resulting  matrix  problem  is: 


M, 


(e)dqj 


+  li1  I  )jj  ([j  =  B'ij'qj 


(e) , 


y 


dt 


-  (e) 


Qi 
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where  is  the  mass  matrix,  is  the  differentiation  matrix  (a  discrete  repre¬ 
sentation  of  V),  is  the  boundary  matrix,  and  is  the  Laplacian  matrix.  The 
discretized  matrices  are: 


(e)  _ 


Mn 


=  Y1u}i 


1=1 


vTn[e) 

uk  u ij 


Mn 

1=1 


Mn 

k=l 


Q+l 

B-f  =  ^2  u>i\Ji\  ipijfi  • 
1=1 


Mn 

J2U1  \Ji\ 


;=i 


where  Q  represents  the  quadrature  used  for  evaluation  of  the  integrals,  which  when 
using  inexact  integration  is  equal  to  N.  Direct  stiffness  summation  (DSS)  was  then 
used  to  construct  the  global  problem: 


Ne 

MU  =  A  Mr 


(e) 

V 


e=l 


LlJ  — 


/\L 


(e) 


e=l 


Djj  = 


A  D 


(e) 

ij 


e=l 


where  Mjj  reduces  to  a  diagonal  matrix  Mj  when  using  inexact  integration. 
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1.  Equation  Set  1 

Applying  these  operators  to  the  Navier-Stokes  Equation  Set  1  (2. 2-2. 3)  yields: 


^  -  ^jMy' DJjUj 

r\  —> 

-yyj-  =  -vTj  My' Duuj  -  cp0Myl DU7Tj  -  gMy'k 

^  =  -ujMy1DIJeJ 


2.  Equation  Set  2 

Similarly  to  Set  1,  the  operators  applied  to  the  Navier-Stokes  Equation  Set  2 
(2. 4-2. 6)  yield: 


dpi 

dt 


~Mj  1  Djj(pu)j 


dui 

~dt 


-ujMj  1Dijuj  -  —Mj  1DijPj  -  gMj  lk 
Pi 


d6j 

~dt 


—vTjMyxDjj6j 
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D.  TEMPORAL  DISCRETIZATION  RUNGE-KUTTA  35 

Integrating  nonlinear  PDEs  forward  in  time  can  lead  to  numerous  problems 
(i.e.  spurious  oscillations,  overshoots,  progressive  smearing,  etc...)  if  the  proper  time 
integration  scheme  is  not  used  [8] .  For  the  proposed  test  cases,  a  strong  stability  pre¬ 
serving  (SSP)  Runge-Kutta  (RK)  method  was  selected  and  implemented  as  outlined 
in  Ruuth  and  Spiteri  [8].  SSP  time  integration  methods  have  strong  nonlinear  stabil¬ 
ity  properties,  which  make  them  optimal,  for  this  thesis,  for  temporal  discretization 
because  of  the  nonlinearities  present  in  the  Euler  and  Navier-Stokes  equations.  The 
particular  method  used  for  time  integration  was  an  explicit  third-order  five-stage  RK 
method  (RK35).  This  method  was  chosen  for  its  large  stability  region  relative  to 
other  explicit  methods  of  equal  order.  The  RK35  method  can  be  represented  by: 


<?„  =  Qin) 

I- 1 

Qi  =  X!  (ai,kQk)  +  1  —  1,2, ...,  s 

k= 0 

Qn  +  1  =  q(s) 

where  I  is  the  stage  (5  for  RK35)  and  the  coefficients  atJ  and  $  are  listed  in  the 
Appendix. 
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III.  APPLIED  COORDINATE 
TRANSFORMS 

A.  EQUATION  SET  1 

1.  Perturbation  Method 

This  section  will  outline  the  expansion  of  the  terms  (application  of  the  pertur¬ 
bation  method)  of  set  1  of  the  Navier-Stokes  Equations,  where  both  n  and  9  will  be 
split/decomposed  into  two  components,  the  mean  values  (tf  and  6)  and  their  associ¬ 
ated  perturbations  (ir'  and  O')  such  that: 


7 r  =  7r(z)  +  7 r'(x,  z,  t ) 


and 


6  =  6(z)  +  9'(x,  z,  t ) 


where  the  mean  values  satisfy  a  hydrostatically  balanced  atmosphere.  After  lineariza¬ 
tion,  the  pressure  tendency  Eq.  (2.2)  becomes: 


d(n  +  7T )  _  ..  R,  ^ 

- - -  +  u  ■  V(7r  +  n‘ )  -I - (7 r  +  7r  )V  •  u  —  0 

at  cv 


which,  after  simplification,  becomes: 


di r  _  _  .  <9tt  i?,  ^ 

— - h  U  ■  V7 r  +  -U7— - 1 - (7T  +  7T  )V  •  U  =  0 

at  oz  cv 


(3.1) 
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Expansion  of  the  terms  of  the  momentum  Eq.  (2.2)  yields: 


du 

~dt 


+  u  ■  Vu  +  cp(9  +  O') V (7 r  +  7 r') 


which  can  be  expanded  to: 


du 

~dt 


~\~  u  ■  V u  -\-  Cp{6  -\-  O' ’) 


/  chr  cbr  \  /  cbr'  cbr'  \ 
\  dx  ’  J  y  dx  ’  dz  ) 


from  which,  using  the  hydrostatic  approximation  =  —-f-§  yields: 


du 

~dt 


-f  m  ■  Vu  "f  cp(6  O' ) 


g  y  f  dn'  dn'\ 

^ ok  +  yih’ih) 


which  can  be  simplified  to: 


du 

~dt 


+  u  ■  Vu  +  cp(0  +  O') 


f  dir'  c^tt'  \ 
\  dx  ’  dz  ) 


-gk-gjk 


and  finally  becomes: 


du 

~dt 


+  u  •  Vu  +  cp(0  +  0')V n' 


a°-k 

9r 


(3.2) 


Expansion  of  the  terms  of  the  thermodynamic  energy  Eq.  (2.3)  yields: 


d(0  +  O') 
dt 


+  u-  V(0  +  6') 


0 
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which  further  leads  to: 


dO'  _  _  80  n 

—  +  u  ■  W  +  w—  =  0 
at  dz 


(3.3) 


2.  Transform 

Using  the  basic  machinery  prescribed  in  Eqs.  (2.13)  -  (2.22),  the  set  of  non¬ 
conservative  Navier-Stokes  Eq.  (3.1  -  3.3)  were  transformed  from  x-z  coordinates 
to  x-cr~  coordinates.  The  first  machinery  applied  to  the  pressure  tendency  Eq.  (3.1) 
was  to  change  the  vector  dot  products  to  transposes  (i.e.  u- V  to  uTV),  which  yielded: 

dir'  _/r,  ,  dir  R, 

— — I-  u 1  X7ir  +  w— — ( - (it  +  7r  )V  ■  u  —  0 

at  oz  cv 

substituting  in  Eq.  (2.13)  and  Eq.  (2.14)  leads  to: 

% + + £<* + v)<  j-vf  (j-aj = o 

applying  the  linear  algebra  identity  in  Eq.  (2.15)  and  the  transformation  for  ( w ) 
yields: 


<9tt' 

~dt 


(Z)T(rl)T{jT)(y)ir'  + 


(  ~dzs  a z-  H 
{  U  dx  H 


+  w 


H-zs\ 

H  ) 


dir  daz 
daz  dz 


+  ^(vr  +  tt')(V)t(Jt)t(J  x)(u)  =  0 
cv 
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which  simplifies  to: 


dir' 

~dt 


(^)T(V)7r'  + 


(  ~dzs  crz  —  H 
{  Udx  H 


+  w 


H-zs\ 

H  J 


dir  /  H  \ 
daz  \H  —  zs) 


+  —  (7T  +  7r')(V)T(-u) 


0 


and  then  to: 


dir'  ~  ~  , 

— +  U.VX 


u 


Or. 


H  \  dzs  dir  „  dir  R,  .  ~  ~ 

—  ]  —  W - 1-  W- - h  —  (IT  +  7T  )V  •  U  =  0 

OOr  Cv 


H  —  ZeJ  dx  da- 


and  then  further  simplified  to: 


dir  ~  ~  , 

— - h  U  ■  V  7T  —  U 

dt 


(Jz-H 

H-z< 


dzs  dir  dir  R,  ~ 

+  - 1 - (vr  +  7r)V-n  =  0 

acr,  ax  acr,  cv 


to  finally  yield: 


dir  ~  ~  .  dir  R ,  .  ~  ~ 

— - hM'V7T  +iy— - 1 - (7T  +  7T  )V  •  «  =  0  (3.4) 

dt  doz  cv 


Applying  the  same  machinery  as  above  to  the  momentum  Eq.  (3.2),  the  dot  products 
were  replaced,  which  yields: 


— 

—  +  uTVu  +  cp(6  +  6')Vir' 


c9-k 

9e 


then  u  and  V  were  replaced  by  their  transforms  (Eq.  (2.13)  and  Eq.  (2.14))  leading  to: 

dTl  _  f 

^  +  (J-^)T(JTV)n  +  cp(9  +  9'){JTV)ir'  =  gjk 
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applying  the  linear  algebra  identity  in  Eq.  (2.15)  yields: 


^  +  (^)T(J-1)T(JT)(V)n  +  cp(9  +  9'){JTvy  =  gjk 


which  simplihes  to: 


^  +  (u)T(V)u  +  cp(9  +  6')(JTvy  =  g^k 


to  finally  yield: 


r)it  ~  ~  fr  —> 

—  +  u  •  Vu  +  cp{9  +  6>')(JT V)7r'  =  gjk 


Applying  the  machinery  to  thermodynamic  energy,  Eq.  (3.3)  yields: 


09’  09  n 

—  +  uTV9  +  w—  =  0 
at  Oz 


substituting  in  Eq.  (2.13)  and  Eq.  (2.14)  leads  to: 


+  (J~1u)T(JTV)d/  +  w 


06  0az 
Oz  0oz 


applying  the  linear  algebra  identity  in  Eq.  (2.15)  and  the  transformation  for  (w) 
yields: 


+  (u)T(J-1)T(JT)(V)d/+  \  -u 


Oz ,  or  —  H  —  zq\  09  Oo 


Ox  H 


H  /  0oz  Oz 
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which  simplifies  to: 


d&_ 

8t 


+  (If)T(V)d'  + 


(  ~dzs  az  —  H 
{  udx  H 


+  w 


H-zs\ 

H  ) 


89  /  II  \ 

8az  \H  —  zs) 


0 


and  simplifies  further,  yielding: 


89'  ~  ~  .  „crz  —  H  8zs  89 

—  +  u.\79'-u- - - 

8t  H  —  zs  8x  8az 


+  w 


89 

daz 


0 


and  then  to: 


89'  ~  ~  .  _az  —  H  8zs  89  89 

_  +  K.W 


0 


to  finally  yield: 


89’  ~  -  _  89  n 

—  +  u-W9  +  w— —  =  0 
at  8az 


(3.6) 


3.  Decomposition 

In  order  to  discretize  the  governing  equations  and  code  them  into  Fortran 
90/95,  the  vector  fields  had  to  be  decomposed  into  scalar  components.  The  decom¬ 
position  of  the  pressure  tendency  Eq.  (3.4)  became: 


8i r' 

'  dir'  dii'' 

dn  R  . 

8u 

8w  1 

~8t  + 

u  a  +  w 

8x  8az 

+  +  (vr  +  7T  ) 

8<jz  cv 

8x 

- 1 

b 

0 
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(3.7) 


The  momentum  Eq.  (3.5)  decomposes  into  two  separate  equations  (u  and  w),  with 
the  first  equation  (u)  taking  the  form: 


du 

dt 


_  du  _  du 

U-  +  W-— 

dx  daz 


+  Cp(9  +  9 ') 


did_ 
dx  + 


(  dzs  az-  H\ 
\  dx  H  —  zs ) 


du' 

daz 


0 


(3.8) 


and  the  representation  for  (w)  taking  the  form: 


dw 


+  w 


dw 

daz 


+  Cp{9  +  6') 


\(  H  ) 

\H  —  zs) 

- 1 

t'* 

b 

<^> 

(3.9) 


The  decomposed  thermodynamic  energy  Eq.  (3.6)  yields: 


d&_ 

~dt 


'  dO'  d9' 

u—  +  w~— 

dx  da? 


+  w 


d6 
da , 


=  0 


(3.10) 


4.  Application  of  the  Galerkin  Statement 

Using  the  Galerkin  machinery  outlined  in  the  previous  chapter,  Eq.  (3.7)  - 
(3.10)  were  discretized  yielding  the  continuous  Galerkin  Set  1  using  x-az  coordinates 
(CGI  x-az): 


Mn  dn'- 

E^  I  Ji  I  A,iA,i  wr 

1=1  UT 

Mn 

E^  \Jl\  A, l 

1=1 


Mn 

E 

k=\ 


I  ( ~  dlpj  i  ,  ~  dlpj  i  ,  _  &lpj  i  _ 


R 


(mn 


- E  + 


\k=l 


dx 


Uj  + 


dipj,i 

da. 


(3.11) 
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Mjv  du  ■ 

l=i 


dt 


mn 

1=1 


-cP(jTMh+0'k)) 


dzs\  &z,j  H 

,  <9^  J  j  H  -  zSyj 


d^i> 
daz  j 


Mn  dw  • 

i=i 


dt 


mn 

\M  A, i 


1=1 


,  ( ~  dtl>j,i  .  d^jtl 

-  lJ>w  \Uk~^rwj +  Wk^rwi 


d  crz 


' mn 

-cP  ( 53^(0*  +  ^) 


» /c=i 


H 


H  Zsj 


— ^  +07T 


dcr. 


0< 


Mn  dd' 

J2ui  \M  = 

1=1  OT 


mn 

\Ji\i>i,i 


i=i 


Mn 


fc=l 


The  existing  CGI  x-z  code  (in  Fortran  90/95)  [1]  was  modified  using  Eqs. 
(3.14)  and  then  used  for  the  test  cases  outlined  in  the  next  chapter. 


(3.12) 

(3.13) 

(3.14) 

(3.11)  - 


B.  EQUATION  SET  2 


1.  Perturbation  Method 

This  section  will  outline  the  expansion  of  the  terms  of  set  2  of  the  Navier- 
Stokes  Equations,  where  (as  seen  in  set  1)  both  p  and  6  will  be  split/decomposed  into 
two  components,  the  reference  values  (p  and  6)  and  their  associated  perturbations  {p1 
and  6').  As  seen  in  set  1: 


p  =  p(z)  +p\x,z,t) 


and 


6  =  6{z)  +  0'(x,  z,  t ) 


where  the  reference  values  again  satisfy  a  hydrostatically  balanced  atmosphere.  In 
order  to  linearize,  the  mass  Eq.  (2.4)  the  product  rule  is  first  applied,  yielding: 


dp  _ 

—  +  u  ■  Vp  +  pV  ■  u 
at 


0 


followed  by  expansion  of  the  terms,  becomes: 

p,)  +  a  •  v(p + p')  +  (p + p'y v  •  u = o 


which,  after  simplification,  yields: 
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(3.15) 


%■  +  u  •  Vp'  +  w^-  +  (p  +  ,/)V  ■  u  =  0 
at  az 


By  setting  p  =  0,  the  momentum  Eq.  (2.6)  becomes: 


1  -» 
—  +  u  •  W  +  -  VP  =  —  gk 
at  p 


and  followed  by  an  expansion  of  the  terms,  yields: 


Ou  1  —  -» 

m+a-va+JpT?)V{p+p,)  =  -3k 


which  leads  to: 


du  _  1  ,  1  d  Pr  - 

a+“'v“+(yM)VF  +  (V7)&i;  =  “9A: 


from  which,  using  the  hydrostatic  approximation  =  —pg  yields: 


du  _  1  pq  7  r 

m+u'Vu  +  WT?)Vp-WT?)k  =  -sk 


which,  after  simplification,  becomes: 


^  +  «-v»+  1  -VP/+  f/9  fc  =  0 

(p  +  p')  (p  +  p7) 
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and  finally: 


du  _  __  1  o' a  T 

—  +  u  ■  Vu  +  — - -VP  =  - —  k 

dt  (p  +  p')  (p  +  p') 


(3.16) 


By  setting  p  =  0,  the  thermodynamic  energy  Eq.  (2.6)  becomes: 


^  +  u  •  W  =  0 

dt 


and  followed  by  an  expansion  of  the  terms,  yields: 


d  (9  +  9') 
dt 


+  u-X7{6  +  6')  =  0 


which  further  leads  to: 


89'  _  _  n 

—  +  u-  V0  +  w—  =  0 
dt  dz 


(3.17) 


2.  Transform 

Using  the  basic  machinery  described  in  Eq.  (2.13)  -  (2.22),  the  set  of  non¬ 
conservative  Navier-Stokes  Eq.  (3.15  -  3.17)  were  transformed  from  x-z  coordinates 
to  x-oy  coordinates.  The  first  machinery  applied  to  the  pressure  tendency  Eq.  (3.15) 
was  to  change  the  vector  dot  products  to  transposes  (i.e.  u- V  to  nTV),  which  yielded: 

+  u'Vp1  +  +  (p  +  p')S7ru  =  0 

dt  dz 
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substituting  in  Eq.  (2.13)  and  Eq.  (2.14)  leads  to: 


dp' 

~dt 


+  {J-l~u)T(JT\7)p'  +  w 


dpdaz 
dz  daz 


+  (p  +  p')(JTV)T  J-^u 


0 


applying  the  linear  algebra  identity  in  Eq.  (2.15)  and  the  transformation  for  (w) 
yields: 


dp' 

~dt 


+  (iff  (j-y^vy  + 


(  „dzscrz-  H 
{  U  dx  H 


+  w 


H  —  zs\ 

H  ) 


dp  daz 
daz  dz 


+  (p  +  p')(V)T(Jr)rJ“1ff  =  0 


which  simplifies  to: 


dp' 

~dt 


+  (tifvp  + 


_ dzs  az  —  H  „  H  —  zs\  dp  (  H  \ 
U~dx  H  +  W  H  )  fcTz  \H-zJ 


+  (p  +  p')(V)TS 


0 


to  finally  yield: 


dp' 

~dt 


+  u  ■  Vp'  +  w 


dp 

daz 


+  (p  +  p') V  •  u  —  0 


(3.18) 


Applying  the  same  machinery  as  above  to  the  momentum  Eq.  (3.16),  the  dot  prod¬ 
ucts  were  replaced,  which  yields: 


du  ->  1 

—  +  M  Vu  +  - - - 

dt  (p  +  pO 


vp'  = 


p'g 


C p  +  p ' 
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then  u  and  V  were  replaced  by  their  transforms  (Eq.  (2.13)  and  Eq.  (2.14))  leading  to: 


^  +  (J~1~u)TJTVu+  1  ,.JTVP'  =  -  P9  k 

dt  (p  +  p')  (p  +  p') 


applying  the  linear  algebra  identity  in  Eq.  (2.15)  yields: 


%  +  ~u(J~l)TJTWu  +  - - - 

dt  y  ’  (p  +  p' 


1  JT\/p'  =  --^—k 


(p  +  p' 


which  simplifies  to: 


84  +  if  VS  +  T^rvp'  =  -+k+k 

dt  (p+p')  (p+p') 


to  finally  yield: 


du  ~  ~  ^  1  tT~  _  p'g  r 

—  +  u  ■  Vu  +  — - -JTVP  =  - —k 

dt  (p  +  p')  (p  +  pf) 


Applying  the  machinery  to  thermodynamic  energy  Eq.  (3.17)  yields: 


d&  ~  d9 

dt  dz 


substituting  in  Eq.  (2.13)  and  Eq.  (2.14)  leads  to: 


d9'  ,  ~  T  T  ~  .  d9  d<Tz 

—  +  (J~lu)T  JTV  9  +  w— — ^  =  0 
dt  dz  duz 


(3.19) 


33 


applying  the  linear  algebra  identity  in  Eq.  (2.15)  and  the  transformation  for  (w) 
yields: 


86' 

~dt 


(  „  8zs  a z  —  H 
{  Udx  H 


+  w 


H-zs\ 

H  ) 


89  8crz 
8az  8z 


0 


which  simplifies  to: 


86' 

~8t 


(tifvd' 


_8zsaz-  H  ~  H  —  zs\  86  (  H  \ 
U~8x  H  +  W  H  )  \H-zJ 


0 


to  finally  yield: 


86'  ~  _ 
— -  +  u  ■  \79  +w 
8t 


86 

8az 


0 


(3.20) 


3.  Decomposition 

Similar  to  CGI,  in  order  to  discretize  the  governing  equations  and  code  them 
into  Fortran  90/95,  they  had  to  be  decomposed  into  scalar  components.  The  decom¬ 
position  of  the  mass  Eq.  (3.18)  became: 


dp'  L  dp'  _  dp' 
~8t  +  U8^z+Wfo- 


+  w 


dp 

8az 


+  (P  +  P') 


8u  8w 
8x  8az 


0 


(3.21) 


The  momentum  Eq.  (3.19)  decomposes  into  two  separate  equations  (u  and  w),  with 
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the  first  equation  (u)  taking  the  form: 


du 

u—  +  w 
dx 


du 

daz 


1 

(P  +  P') 


dP'  (dzs  az  -  H\  dP' 
dx  +  \dx  H  —  zs  )  daz 


0 


(3.22) 


and  the  representation  for  (w)  taking  the  form: 


p'g 

(P  +  P') 

The  decomposed  thermodynamic  energy  Eq.  (3.20)  yields: 


dw 


_  dw  _ dw 

u—  +  w~— 

dx  da7 


+ 


(P  +  P') 


H 


H  —  zj  da 


dP' 


(3.23) 


dP_ 

dt 


'  d0'  d9 ' 

u—  +  w— 
dx  da z 


+  w 


d6 

daz 


=  0 


(3.24) 


4.  Application  of  the  Galerkin  Statement 

Using  the  Galerkin  machinery  outlined  in  the  previous  chapter,  Eqs.  (3.21)  - 
(3.24)  were  discretized  yielding  the  continuous  Galerkin  Set  2  using  x-az  coordinates 
(CG2  x-az): 


mjv  Qp', 

i=i  at 


mn 

\Ji\  A,i 

1=1 


Mn 

A,l  I  A 


k=l 


djjj.i 

dx 


+ 


'  dZe 


a 


z,3 


dx  )  .  H 


H\  dA,i 

da7 


's,j 


+Wk 


H 


H  -  zSjj 


dipj,i  ,  _  dipj,i  -  N 


(mn  \  / 

-  ( J2MPk  +  p'k)  j  y 


dip 


dx 


uj  + 


dip 


3,1 


da , 


Wj 


(3.25) 
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The  existing  CG2  x-z  code  (in  Fortran  90/95)  was  modified  using  Eq.  (3.25)  -  (3.28) 
and  then  used  for  the  test  cases  outlined  in  the  next  chapter. 
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IV. 


TEST  CASES 


A.  OVERVIEW 

Three  specific  cases  where  chosen  to  evaluate  the  a-coordinate  transforma¬ 
tions:  the  rising  thermal  bubble,  the  linear  hydrostatic  mountain,  and  the  linear  non¬ 
hydrostatic  mountain.  The  rising  thermal  bubble,  case  1,  was  chosen  as  a  benchmark 
in  order  to  insure  that  in  the  absence  of  terrain  the  model  dynamics  still  functioned 
properly  (i.e.  the  cr-coordinate  transformed  equations  reduce  to  the  original  govern¬ 
ing  equations)  yielding  the  same  results  as  the  unmodified  source  code  [1].  The  linear 
hydrostatic  mountain,  case  2,  and  the  linear  non- hydrostatic  mountain,  case  3,  were 
chosen  to  evaluate  the  advantages  and  disadvantages  of  the  x-cr  governing  equation 
in  relation  to  their  x-z  formulation  in  both  hydrostatic  and  non-hydrostatic  environ¬ 
ments  using  the  same  terrain.  This  chapter  will  outline  the  test  case  assumptions, 
additional  information  derived  for  the  cr-coordinates,  and  previous  results  that  will 
be  used  for  comparison. 

B.  CASE  1:  RISING  THERMAL  BUBBLE 

This  test  case  represents  a  highly  nonlinear  non- hydrostatic  flow  problem,  by 
demonstrating  the  evolution  of  an  initially  at  rest  warm  air  bubble,  relative  to  the 
surrounding  environment,  in  a  hydrostatic  constant  potential  temperature  environ¬ 
ment  as  resolved  by  the  Navier-Stokes  equations.  As  the  warm  air  bubble  rises,  it  will 
deform,  taking  the  shape  of  a  mushroom  cloud,  clue  to  the  shearing  motion  caused  by 
the  resulting  differential  vertical  velocity  field.  The  initial  distribution  of  potential 
temperature  perturbations  in  order  to  drive  vertical  motion  is  defined  by: 


O' 


0 

8c 

2 


1  +  COS 


r  >  r c, 
r  <  rc, 
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where  9C  =  0.5°C,  ttc  is  the  trigonometric  constant,  r  =  \J {x  —  xc )2  +  (z  —  zc )2  with 
the  constants:  6  =  300  K,  rc  =  250  m,  and  (xc,  zc)  =  (500,350)  m.  [1]  The  domain 
of  interest  is  (x,z)  G  [0, 1000]2  m  integrated  over  the  time  interval,  t  G  [0,700]  s. 
Additionally,  the  boundary  conditions  for  all  boundaries  are  no-flux. 

The  utility  of  running  this  case  is  to  insure  that  the  a  transformed  governing 
equations  yield  the  same  results  as  the  x-z  set  of  the  governing  equations  when  the 
terrain  is  flat: 


Zsurf  0 


where  zsurf  is  the  height  of  the  surface.  The  surface  heights  are  vital  in  deriving  the 
cr-coordinates  (see  equation  2.16).  The  slope  of  the  surface  must  also  be  calculated 
for  the  transformation,  which  in  this  case  leads  to: 

d  Zsurf 

—  =  0 

With  a  flat  surface  in  combination  with  a  slope  of  0,  the  transformation  of  the  coor¬ 
dinate  system  (x-za)  reduces  to  x-z,  making  both  the  modified  and  unmodified  codes 
yield  identical  results. 

The  unmodified  source  code  was  run  in  order  to  establish  a  baseline  to  evaluate 
the  transformed  source  code.  Three  resolutions  were  used:  20  m  (2601  grid  points),  10 
m  (10201  grid  points),  and  5  m  (40401  grid  points).  The  resulting  data  for  potential 
temperature  perturbations  was  then  plotted  for  each  resolution  after  700  s  of  model 
integration  (see  Figures  2  and  3).  As  seen  in  Figures  2  and  3,  the  warm  air  bubble 
did  deform  into  a  mushroom  type  formation,  while  maintaining  symmetry. 
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Figure  2.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  pertur¬ 
bations  using  unmodified  CGI  source  codes  [1]  after  700  s  for  resolutions:  (a)  20,  (b) 
10,  and  (c)  5  m.  All  cases  were  run  using  10th  order  polynomials,  with  contours  from 
-0.05  to  0.525  K  with  an  interval  of  0.025  K. 
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Figure  3.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  pertur¬ 
bations  using  unmodified  CG2  source  codes  [1]  after  700  s  for  resolutions:  (a)  20,  (b) 
10,  and  (c)  5  m.  All  cases  were  run  using  10th  order  polynomials,  with  contours  from 
-0.05  to  0.525  K  with  an  interval  of  0.025  K. 
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C.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN 


This  test  case  was  chosen  in  order  to  evaluate  the  model  performance  when 
a  simple  terrain  feature  is  introduced  in  a  linear  hydrostatic  flow  environment  with 
steady  inflow  and  outflow  lateral  boundary  conditions.  Over  time,  a  steady-state 
mountain  wave  over  a  single  peak  should  form  if  the  model  dynamics  are  resolving 
the  scenario  accurately.  Initially,  the  constant  mean  horizontal  velocity  is  set  to  u  — 
20  m/s,  the  mean  vertical  velocity  is  set  to  w  =  0  m/s,  and  the  atmosphere  is  set 
to  isothermal  with  a  constant  mean  temperature  of  T  —  250  K.  [1]  Since  this  case 
is  isothermal,  the  buoyancy  frequency  or  Brunt- Vaisala  frequency  J\f  2  =  g4;(ln9) 
reduces  to  A f  —  — /=,  which  simplifies  the  Exner  pressure  to: 


9 

7 r  =  e  cpt 


The  domain  of  interest  is  (x,z)  G  [0,  240,000]  x  [0,  30,000]  m  integrated  over  the  time 
interval,  t  G  [0, 10]  h.  For  this  test  case,  the  terrain  (the  versiera  cli  Agnesi  mountain 
profile)  is  represented  by: 


he 

Zsurf =  7  7  72V 

0  +  Is?) ) 

where  zsurf  is  the  height  of  the  surface,  and  hc,  xc,  and  ac  are  constants  ( hc  =  1  m,  xc 
=  120,000  m,  and  ac  =  10,000  m).  Additionally,  the  bottom  boundary  conditions  are 
no-flux,  while  the  top  and  lateral  boundaries  use  a  non-reflecting  boundary  condition 
[1],  For  further  information  on  such  boundary  conditions,  see  the  papers  by  Durran 
and  Klemp  [9],  Giraldo  and  Restelli  [1],  and  Dea,  Giraldo,  and  Neta  [10].  However 
non-reflecting  boundary  conditions  are  outside  the  scope  of  this  thesis  and  will  not  be 
mentioned  further.  Unlike  the  previous  case,  the  surface  heights  do  vary  with  x  and 
are  even  more  vital  in  deriving  the  cr-coordinates.  The  surface  plot  can  be  seen  in 
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Figure  4. a,  on  which  it  can  be  seen  (noticing  the  z-axis  scale)  that  the  peak  is  fairly 
low  and  has  smooth  geometry.  Since  the  surface  height  does  vary  with  x,  the  slope  of 
the  surface  must  also  be  calculated  for  the  transformation,  which  in  this  case  leads  to: 


Zsurf  he  (  1 


X  —  XC 


ac 


-l 


Zsurf  he 


(ac)2  +  (x2  —  2  (xc)(x)  +  ( xc )2) 


ac 


-l 


zSurf  =  (hc)(acY  (ac)2  +  x2  —  2 (xc)(x)  +  (xc) 


i-i 


dzgurf  _  (—l)(hc)(ac)2(2x  -  2 (xc)) 

dx  [(ac)2  +  x2  —  2  (xc)(x)  +  (xc)2]2 

dzsurf  (—2  )(hc)(ac)2(x  —  xc) 

dx  [(ac)2  +  (x  —  xc)2}2 

dz 

where  the  surface  slope  (  guJf )  plot  can  be  seen  in  Figure  4.b,  on  which  it  can  be 
seen  (noticing  the  z-axis  scale)  that  the  slope  of  the  peak  is  fairly  flat  and  again  has 
smooth  geometry. 

The  unmodified  source  code  was  run  (colored  lines)  and  compared  with  the 
analytic  solutions  (black  lines),  generated  for  the  unmodified  source  code  and  de¬ 
veloped  by  Giraldo  and  Restelli  [1],  in  order  to  establish  a  baseline  to  evaluate  the 
transformed  source  code.  One  resolution  was  used:  1200  m  (in  x)  and  240  m  (in  z) 
(20301  grid  points).  The  resulting  data  for  vertical  and  horizontal  velocity  was  then 
plotted  for  the  resolution  after  10  h  (see  Figure  5).  As  seen  in  Figure  5,  a  steady-state 
mountain  wave  did  form  in  the  appropriate  region  and  will  be  compared  against  the 
transformed  governing  equation  sets. 
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Figure  5.  Case  2:  Linear  Hydrostatic  Mountain.  Resulting  horizontal  velocity  (a)  and 
vertical  velocity  (b)  using  unmodified  (i)  CGI  and  (ii)  CG2  source  codes  [1]  after  10 
h  for  the  resolution  of  1200  m  (in  x)  and  240  m  (in  z)  (colored  lines)  plotted  with  the 
analytic  solution  (black  lines).  All  cases  were  run  using  10th  order  polynomials,  with 
contours  from  -0.025  to  0.025  ms”1  with  an  interval  of  0.005  ms”1,  (a),  and  -0.005  to 
0.005  ms”1  with  an  interval  of  0.0005  ms”1,  (b). 
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D.  CASE  3:  LINEAR  NON-HYDROSTATIC  MOUNTAIN 

The  final  test  case  was  chosen  in  order  to  evaluate  the  model  performance  when 
a  simple  terrain  feature  is  introduced  in  a  linear  non-hydrostatic  flow  environment 
with  steady  inflow  and  outflow  lateral  boundary  conditions.  Over  time,  a  steady-state 
mountain  wave  over  a  single  peak  should  form  if  the  model  dynamics  are  resolving  the 
scenario  accurately.  Initially,  the  constant  mean  horizontal  velocity  is  set  to  u  —  10 
m/s,  the  mean  vertical  velocity  is  set  to  w  =  0  m/s,  and  the  atmosphere  is  uniformly 
stratified  with  a  Brunt-Vaisala  frequency  of  Af  =  0.01  /s.  This  test  case  uses  the 
same  terrain  as  case  2  (see  Figure  4).  The  constants  used  for  this  case  were  hc  =  1 
m,  xc  =  72,000  m,  ac  =  1,000  m,  and  6$  =  280K.  The  domain  of  interest  is  (x,z)  e  [0, 
144,000]  x  [0,  30,000]  m  integrated  over  the  time  interval,  t  G  [0,  5]  h.  Like  case  2, 
the  bottom  boundary  conditions  are  no-flux  and  the  top  and  lateral  boundaries  are 
non- reflecting.  [1] 

The  unmodified  source  code  was  again  run  (colored  lines)  and  compared  with 
the  analytic  solution  [1]  (black  lines)  in  order  to  establish  a  baseline  to  evaluate  the 
transformed  source  code.  One  resolution  was  used:  360  m  (in  x)  and  300  m  (in  z) 
(40501  grid  points).  The  resulting  data  for  vertical  and  horizontal  velocity  was  then 
plotted  after  5  h  (see  Figure  6).  As  seen  in  Figure  6,  a  steady-state  mountain  wave 
did  form  in  the  appropriate  region  and  will  be  compared  against  the  transformed 
governing  equation  sets. 
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Figure  6.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Resulting  horizontal  velocity 
(a)  and  vertical  velocity  (b)  using  unmodified  (i)  CGI  and  (ii)  CG2  source  codes  [1] 
after  5  h  for  the  resolution  of  360  m  (in  x)  and  300  m  (in  z)  (colored  lines)  plotted  with 
the  analytic  solution  (black  lines).  All  cases  were  run  using  10th  order  polynomials, 
with  contours  from  -0.025  to  0.025  ms'1  with  an  interval  of  0.005  ms”1,  (a),  and 
-0.005  to  0.005  ms”1  with  an  interval  of  0.0005  ms”1,  (b). 
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V. 


RESULTS 


A.  CASE  1:  RISING  THERMAL  BUBBLE 

1.  Accuracy  and  Comparison 

The  modified  source  code  (containing  the  cr-coordinate  transformed  governing 
equations)  was  run  using  the  same  three  resolutions  that  were  used  for  the  unmodified 
code:  20  m  (2601  grid  points),  10  m  (10201  grid  points),  and  5  m  (40401  grid  points). 
The  resulting  data  for  potential  temperature  was  then  plotted  for  each  resolution  after 
integrating  forward  700  s  (see  Figure  7).  As  seen  in  Figure  7,  the  warm  air  bubble 
did  deform  into  a  mushroom  type  formation,  while  maintaining  symmetry,  appearing 
to  duplicate  the  results  seen  by  the  x-z  code  (see  figure  2). 

To  further  evaluate  the  similarities  in  the  x-z  and  x-az  model  solutions  nine 
variables  (using  5  m  resolution)  were  compared:  the  maximum  Exner  pressure  pertur¬ 
bation  ( n'max  (unitless)),  the  minimum  Exner  pressure  perturbation  (n'min  (unitless)), 
the  maximum  horizontal  wind  velocity  ( umax  (ms-1)),  the  minimum  horizontal  wind 
velocity  (umm  (ms-1)),  the  maximum  vertical  wind  velocity  (wmax  (ms-1)),  the  mini¬ 
mum  vertical  wind  velocity  ( wmin  (ms-1)),  the  maximum  potential  temperature  per¬ 
turbation  (0'max  (K)),  the  minimum  potential  temperature  perturbation  (9'min  (K)), 
and  CPU  time  (s).  These  values  were  compare  to  better  observe  the  relative  maxi¬ 
mum  and  minimum  values  occurring  for  each  model  run.  In  order  to  have  consistent 
CPU  constraints  both  x-z  and  x-a2  for  set  1  code  were  run  in  parallel  on  the  same 
machine  (NPS  Math  Department  32  Processor  Apple  Cluster:  Riemann)  and  started 
at  the  same  time,  and  then  repeated  for  set  2.  The  resulting  values  can  be  seen  in  Ta¬ 
ble  I.  All  four  sets  of  code  converge  to  nearly  identical  solutions.  The  only  difference 
are  the  following:  minor  differences  in  the  i Trmax  for  set  1  and  set  2;  and  differences  in 
0'max  and  9'min  between  x-z  coordinates  and  x-az  coordinates  for  both  sets.  The  major 
difference  between  the  models  is  the  CPU  time  (in  seconds),  which  was  consistently 
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slower  for  the  x-az  (  28.24  %  increase  in  time  for  set  1  and  15.58  %  increase  in  time 
for  set  2)  relative  to  their  x-z  counterparts.  This  difference  is  caused  by  the  additional 
x-az  data  being  passed  between  program  functions  and  additional  calculations  being 
performed  (  0(5  *  [NT]  *  [NE]  *  [ MN ]2  *  [8  *  (IV  +  1)  +  26])),  where  NT  is  the  number 
of  time  steps,  NE  is  the  number  of  elements,  and  MN  is  the  number  of  quadrature 
points  or  order  N  + 1).  For  the  5  m  resolution  case,  there  were  5.2101  x  1019  additional 
operations.  The  horizontal  velocities  umax  and  umin  also  indicate  that  symmetry  is 
maintained  for  all  four  model  sets. 

To  further  evaluate  the  performance  of  the  four  model  runs,  the  Root  Mean 
Squared  Error  (RMSE)  was  determined  for  each  run  in  relation  to  the  analytic  solu¬ 
tion  (see  Table  II)  for  the  four  variables:  i r,  u,  ui,  and  6.  The  RMSE  was  calculated 
using: 


\RMS 


\ 


1 

_  ^  ^  ^ q  numerical  _  g  analytic ^ 

Np  i= i 


where  q  represents  the  given  state  variable  vector,  Np  is  the  number  of  points,  and  the 
analytic  solution  for  this  case  was  defined  by  the  unmodified  source  code  numerical 


Table  I.  Case  1:  Rising  Thermal  Bubble.  Comparison  of  modified  and  unmodified 


CGI  and  CG2  using  5  m  resolution  after  700s. 


Model 

CGI  x-z 

CGI  x-crz 

CG2  x-z 

CG2  x-az 

7 t' 

11  max 

0.9364x  10-5 

0.9364X10-5 

0.9365  xlO-5 

0.9365x  10-5 

7T/  • 
mm 

-0.1195X10-4 

-0.1195X10-4 

-0.1195x  10-4 

-0.1195x  10-4 

UmaX  (mS_1) 

0.2079X101 

0.2079x  101 

0.2079X101 

0.2079X101 

Umin  (niS-1) 

-0.2079X101 

-0.2079X101 

-0.2079X101 

-0.2079X101 

Wmax  (lUS-1) 

0.2536X101 

0.2536x  101 

0.2536X101 

0.2536X101 

Wmin  (lUS-1) 

-0.1912X101 

-0.1912X101 

-0.1912X101 

-0.1912X101 

O'max  (K) 

0.5713x10° 

0.5715x10° 

0.5713x10° 

0.5715x10° 

O'min  (K) 

-0.9736X10-1 

-0.9735xl0-1 

-0.9736x  10-1 

-0.9735x  10-1 

CPU  time  (s) 

0.1027xl05 

0.1317xl05 

0. 1322x10s 

0.1528x10s 
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Table  II.  Case  1:  Rising  Thermal  Bubble  RMSE.  Root- mean-squared  errors  for  the 
four  primary  state  variables,  for  the  modified  codes  using  the  unmodified  code  so¬ 
lutions  as  the  analytic  solutions,  after  700  s  using  5  m  resolution  and  10th  order 
polynomials. 


Model 

IT 

u  (ms  *) 

w  (ms  x) 

*(K) 

CGI  x-az 
CG2  x-az 

2.5931xl0“15 

2.2034xl0”16 

9.2434xl(T10 
5.2908x  10-12 

1.6068xl0-9 

8.5529X10"12 

1.8639x  10-9 
9.6233x  10-12 

solutions  in  order  to  better  compare  x-z  solution  to  the  x-az  solution.  Both  sets  1 
and  2  exhibited  a  smaller  error  for  all  four  variables  of  interest,  with  set  2  having  the 
lower  RMSE  values  overall. 

In  addition  to  comparing  extremes  of  the  selected  fields,  a  one- dimensional 
vertical  temperature  profile  along  the  centerline  (x  =  500  nr)  was  plotted  for  each 
resolution  to  better  discern  the  differences  between  the  various  models  (see  Figure 
9).  As  the  resolution  increased  from  20  m  to  5  m,  there  is  a  noticeable  increase  in 
the  sharpness  of  the  temperature  distributions.  As  expected,  for  all  three  resolutions, 
there  is  no  significant  difference  between  the  temperatures  along  the  centerline  for 
each  model  at  that  specific  resolution. 

2.  Conclusions 

All  four  model  runs  for  the  rising  thermal  bubble  did  deform  into  a  mushroom 
type  formation,  while  maintaining  symmetry  (as  indicated  by  the  distribution  of  hor¬ 
izontal  velocities).  Additionally,  the  models  converged  to  nearly  identical  solutions, 
as  demonstrated  by  the  min  and  max  values  and  the  RMSE.  Overall,  for  both  sets  1 
and  2  of  the  Navier-Stokes  equations,  the  x-crz  coordinates  do  reduce  to  x-z,  showing 
that  in  the  absence  of  terrain  that  the  model  dynamic  still  function  properly,  but  did 
show  a  noticeable  increase  in  computational  expense.  These  results  validate  the  first 
stage  of  the  code  evaluation. 
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Figure  7.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  pertur¬ 
bations  using  modified  CGI  source  codes  [1]  after  700  s  for  resolutions:  (a)  20,  (b) 
10,  and  (c)  5  m.  All  cases  were  run  using  10th  order  polynomials,  with  contours  from 
-0.05  to  0.525  K  with  an  interval  of  0.025  K. 
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1000 


Figure  8.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  pertur¬ 
bations  using  modified  CG2  source  codes  [1]  after  700  s  for  resolutions:  (a)  20,  (b) 
10,  and  (c)  5  m.  All  cases  were  run  using  10th  order  polynomials,  with  contours  from 
-0.05  to  0.525  K  with  an  interval  of  0.025  K. 
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Figure  9.  Case  1:  Rising  Thermal  Bubble.  Resulting  potential  temperature  perturba¬ 
tions  for  all  four  models  along  the  vertical  axis  (x  =  500  m)  after  700  s  for  resolutions: 
(a)  20,  (b)  10,  and  (c)  5  m.  All  cases  were  run  using  10th  order  polynomials. 
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B.  CASE  2:  LINEAR  HYDROSTATIC  MOUNTAIN 


1.  Accuracy  and  Comparison 

The  modified  source  code  was  run  and  the  numerical  solution  after  10  hours 
(colored  lines)  was  plotted  with  and  compared  to  the  analytic  solution  [1]  (black  lines) 
(see  Figure  10).  Only  one  resolution  was  used:  1200  m  (in  x)  and  240  m  (in  z)  (20301 
grid  points).  The  resulting  data  for  the  vertical  and  horizontal  velocity  was  then  used 
for  the  comparison  (see  Table  111).  As  seen  in  Figure  10,  a  steady-state  mountain 
wave  did  form  in  the  appropriate  region  closely  modeling  the  analytic  solution. 

Visually  comparing  Figure  5  to  Figure  10,  indicated  fairly  identical  solutions, 
and  to  further  evaluate  the  similarities  in  the  x-z  and  x-oy  model  solutions  nine 
variables  were  again  compared:  7r'max,  nrmin,  umaxi  umin,  wmax,  wmin,  9'max,  9'min,  and 
CPU  time  (in  seconds).  Similar  to  case  1,  both  x-z  and  x-oy  for  set  1  code  were  run 
in  parallel  and  started  at  the  same  time,  and  then  repeated  for  set  2.  The  resulting 
values  can  be  seen  in  Table  111.  All  four  sets  of  code  converge  to  similar  solutions,  with 
only  minor  differences  overall,  and  the  most  predominate  differences  in  the  wmax  and 
Wmin  between  x-z  coordinates  and  x-oy  coordinates,  where  x-z  had  greater  magnitude 
values  for  wmin  and  x-az  had  greater  magnitude  values  for  wmax.  As  seen  in  case  1, 
another  major  difference  between  the  models  is  the  CPU  time,  which  was  consistently 
slower  for  the  x-oy  (  30.03  %  increase  in  time  for  set  1  and  15.67  %  increase  in  time 
for  set  2)  relative  to  their  x-z  counterparts.  This  difference  is  again  being  caused 
by  the  additional  x-oy  data  being  passed  between  program  functions  and  additional 
calculations  being  performed  (  0(5  *  [jVT]  *  [Ne]  *  [Mat]2  *  [8  *  (N  +  1)  +  26])).  For 
this  case,  there  were  1.8793xl015  additional  operations. 

To  further  evaluate  the  performance  of  the  four  model  runs,  the  R.MSE  was 
determined  for  each  run  in  relation  to  the  analytic  solution  (see  Table  IV)  for  the  four 
variables:  tt,  u,  w,  and  9.  The  analytic  solution  was  interpolated  to  the  model  grid 
using  a  cubic  approximation.  Once  the  domain  of  interest  was  defined  and  the  nu- 
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Table  III.  Case  2:  Linear  Hydrostatic  Mountain.  Comparison  of  modified  and  un¬ 
modified  CGI  and  CG2  using  1200  m  (in  x)  and  240  m  (in  z)  resolution  after  10 
hours. _ 


Model 

CGI  x-z 

CGI  x-oy 

CG2  x-z 

CG2  x-oy 

7 t' 

11  max 

0.1792x  10-5 

0.1786X10-5 

0.1781xl0-5 

0.1785x  10-5 

7T/  • 
mm 

-0.1859X10-5 

-0.1859X10-5 

-0.1840x  10-5 

-0.1846x  10-5 

UmaX  (mS_1) 

0.4072x  10-1 

0.4070X10-1 

0.4034xl0-1 

0.4044x  10-1 

Umin  (ms-1) 

-0.3506X10-1 

-0.3484xl0-1 

-0.3474x  10-1 

-0.3473x  10-1 

Wmax  (niS-1) 

0.4229  xlO-2 

0.4656X10-2 

0.4245X10-2 

0.4686x  10-2 

Wmin  (niS-1) 

-0.5199X10-2 

-0.4279X10-2 

-0.5214x  10-2 

-0.4313x  10-2 

0'max  (K) 

0.2400x  10-1 

0.2386xl0-1 

0.2380xl0-1 

0.2382x  10-1 

0'min  (K) 

-0.3184xl0-1 

-0.3124xl0-1 

-0.3157xl0-1 

-0.3115x  10-1 

CPU  time  (s) 

0.4519X104 

0.5876x  104 

0.6463xl04 

0.7476xl04 

merical  and  analytic  solutions  were  on  matching  grids  a  bootstrap  (random  sampling) 
method  was  used  to  calculate  the  95  %  confidence  interval  (Cl)  [11].  The  95%  Cl 
was  needed  in  order  to  determine  if  the  differences  observed  in  the  RMSE  values  were 
indeed  significant.  The  bootstrap  method  built  the  95%  Cl  by  creating  a  domain  (of 
equal  size  to  the  domain  of  interest)  that  was  populated  with  random  samples  from 
the  original  numerical  and  analytic  solution  pairs,  and  then  the  RMSE  for  the  new 
domain  was  calculated  and  stored.  This  process  was  iterated  10,000  times,  storing 
the  RMSE  from  each  iteration.  The  derived  RMSE  values  were  then  sorted  and  the 
values  at  2.5%  and  97.5%  of  the  distribution  were  taken  and  compared  to  the  original 
RMSE  in  order  to  establish  the  95%  CL  The  resulting  confidence  intervals  indicate 
that  the  differences  between  the  state  variable  RMS  errors  are  significant,  since  the 
ranges  do  not  overlap.  For  set  1,  x-z  coordinates  exhibited  a  significantly  smaller 
error  for  all  four  variables  of  interest.  For  set  2,  x-oy  had  lower  RMSE  values  for  n 
and  u,  but  larger  RMSE  values  for  w  and  6.  Taking  a  closer  look  at  the  w  RMSE 
with  respect  to  wmax  shows  that  CGI  x-z  has  a  1.32%  error  as  compared  to  4.24% 
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Table  IV.  Case  2:  Linear  Hydrostatic  Mountain  RMSE.  Root-mean-squared  errors  for 
the  four  primary  state  variables,  for  both  the  modified  and  unmodified  codes,  after 
10  h  using  1200  m  (in  x)  and  240  m  (in  z)  resolution  and  10th  order  polynomials. 


Model 

7T 

95%  Conhdence  Interval 

CGI  x-z 

1.2263x  10-7 

+  8.9354xlO-10 

-  8.9757xlO-10 

CGI  x-crj 

1.2898x  10-7 

+  9.5565xlO-10 

-  9.9149xl0-1° 

CG2  x-z 

1.2121x  10-7 

+  9.1186xlO-10 

-  9.1327x  10-1° 

CG2  x-az 

1.2200x  10-7 

+  9.0255xl0-1° 

-  9.1227x  10-1° 

Model 

u  (ms-1) 

95%  Conhdence  Interval 

CGI  x-z 

2.2561  xlO-3 

+  1.7930X10-5 

-  1.8157x  10-5 

CGI  x-cr2 

2.4089x  10-3 

+  1.9437xl0-5 

-  1.9878x  10-5 

CG2  x-z 

2.2481x  10-3 

+  1.8488X10-5 

-  1.8337xl0-5 

CG2  x-cr2 

2.2808x  10-3 

+  1.8327xl0-5 

-  1.8985x  10-5 

Model 

w  (ms-1) 

95%  Conhdence  Interval 

CGI  x-z 

6.8468x  10-5 

+  6.6744xl0-7 

-  6.6389x  10-7 

CGI  x-cr2 

1.9757xl0-4 

+  2.1074xl0-6 

-  2.1468x  10-6 

CG2  x-z 

5.1424x  10-5 

+  5.0225X10-7 

-  5.0915x  10-7 

CG2  x-cr2 

1.9311x  10-4 

+  2.0131X10-6 

-  2.0667x  10-6 

Model 

d(K) 

95%  Conhdence  Interval 

CGI  x-z 

1.6565x  10-3 

+  1.3026X10-5 

-  1.3132x  10-5 

CGI  x-cr2 

1.8506x  10-3 

+1.4243x  10-5 

-  1.4748x  10-5 

CG2  x-z 

1.6360x  10-3 

+  1.4812X10-5 

-  1.4410x  10-5 

CG2  x-cr2 

1.7208x  10-3 

+  1.4790X10-5 

-  1.5027x  10-5 

for  CGI  x-a2  and  that  CG2  x-z  has  a  0.99%  error  as  compared  to  4.12%  for  CG2 
x-ox,  which  is  substantial  for  both  sets. 

In  order  to  verify  that  a  steady-state  solution  for  the  linear  hydrostatic  moun¬ 
tain  case  was  achieved,  the  momentum  flux  was  derived  using: 


OO 

m(z)  =  J  p(z)u(x,  z)w(x,  z)dx  (5.1) 

—  OO 


where  p(z)  is  the  reference  density  as  a  function  of  height  (z).  [12]  Additionally,  the 
analytic  hydrostatic  momentum  flux  (from  linear  theory)  is  given  by: 
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(5.2) 


TTc 

4 


psusNh?c 


where  mH  denotes  the  analytic  hydrostatic  momentum  flux,  ps  is  the  reference  density 
at  the  surface,  us  is  the  horizontal  velocity  values  at  the  surface,  A f  is  the  Brunt- 
Vaisala  frequency,  and  hc  is  the  height  of  the  mountain.  [1]  The  momentum  flux  was 
then  normalized  by  m(z)/mH(z),  and  will  be  the  value  discussed  in  this  thesis.  In 
Figure  11,  the  normalized  momentum  flux  for  all  four  model  runs,  using  1200  m  (in 
x)  and  240  m  (in  z)  resolution,  was  plotted  for  2  h,  4  h,  6  h,  8  h,  and  10  h.  As  seen  in 
Figure  11,  all  four  model  runs  converge  to  a  steady  state  solution  after  10  h.  Of  the 
four  model  simulations,  the  CGI  x-z  and  CG2  x-z  models  yielded  results  that  were 
far  better  than  the  x-<r2  models,  with  values  between  0.95  and  1.01  versus  between 
0.65  and  1.15  as  seen  for  the  x-az  runs. 

2.  Conclusions 

All  four  model  runs  for  the  linear  hydrostatic  mountain  case  did  develop  a 
steady-state  mountain  wave  over  a  single  peak  indicating  that  the  model  dynamics 
are  resolving  the  scenario  accurately.  Additionally,  the  models  converged  to  nearly 
identical  solutions.  Overall,  for  both  set  1  and  set  2  of  the  Navier-Stokes  equations, 
the  x-az  coordinates  performed  slightly  worse  than  their  x-z  counterparts  (reflected 
in  the  RMSE,  the  normalized  momentum  flux,  and  most  pronounced  in  w  with  an 
approximately  four  times  larger  RMSE  with  respect  to  wmax  )  and  as  seen  in  the 
previous  case  did  show  a  noticeable  increase  in  computational  expense.  These  results 
appear  to  indicate  that  when  the  model  dynamics  are  resolved  using  purely  explicit 
time  integration  methods  that  there  is  no  clear  advantage  in  using  x-<r2  over  x-z  and 
will  be  further  explored  into  the  next  case. 
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Figure  10.  Case  2:  Linear  Hydrostatic  Mountain.  Resulting  horizontal  velocity  (a) 
and  vertical  velocity  (b)  using  modified  (i)  CGI  and  (ii)  CG2  source  codes  [1]  after  10 
h  for  the  resolution  of  1200  m  (in  x)  and  240  m  (in  z)  (colored  lines)  plotted  with  the 
analytic  solution  (black  lines).  All  cases  were  run  using  10th  order  polynomials,  with 
contours  from  -0.025  to  0.025  ms-1  with  an  interval  of  0.005  ms-1,  (a),  and  -0.005  to 
0.005  ms-1  with  an  interval  of  0.0005  ms-1,  (b). 
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Figure  11.  Case  2:  Linear  Hydrostatic  Mountain, 
the  resolution  of  1200  m  (in  x)  and  240  m  (in  z),  at 
for  the  four  model  runs:  (i)  CGI  x-z,  (ii)  CGI  x-az, 


Normalized  momentum  flux  for 
times  2  h,  4  h,  6  h,  8  h,  and  10  h 
(iii)  CG2  x-z,  and  (iv)  CG2  x-az. 


C.  CASE  3:  LINEAR  NON-HYDROSTATIC  MOUNTAIN 


1.  Accuracy  and  Comparison 

For  the  final  case,  both  CGI  x-oy  and  CG2  x-oy  were  run  and  the  numerical 
solution  after  5  hours  (colored  lines)  was  plotted  and  compared  to  the  analytic  solution 
[1]  (black  lines)  (see  Figure  12).  Similar  to  case  2,  only  one  resolution  was  used:  360 
m  (in  x)  and  300  m  (in  z)  (40501  grid  points).  The  resulting  data  for  the  vertical 
and  horizontal  velocity  was  then  used  for  the  comparison  (see  Table  V).  As  seen  in 
Figure  12,  a  steady-state  mountain  wave  did  form  in  the  appropriate  region,  there 
were  some  apparent  contrasts  with  the  analytic  solution  for  CGI  x-oy  and  CG2  x-oy. 

Visually  comparing  Figure  12.i.a,b  to  Figure  12.ii.a,b  indicated  fairly  similar 
solutions  for  the  horizontal  velocities  (a).  Comparing  the  vertical  velocities  (b)  re¬ 
vealed  a  similar  pattern  in  the  placement  of  the  steady  state  waves  between  x-z  and 
x-oy,  but  the  x-oy  coordinates  appear  to  have  significantly  stronger  oscillations  than 
the  x-z  coordinates.  To  further  evaluate  the  x-z  and  x-oy  model  solutions  nine  vari¬ 
ables  were  again  compared:  7r^aa.,  n'min,  umax,  umin,  wmax,  wmin,  9'max ,  Q'min,  and  CPU 
time  (in  seconds).  Similar  to  case  1  and  case  2,  both  x-z  and  x-oy  for  set  1  code 
were  run  in  parallel  and  started  at  the  same  time,  and  then  repeated  for  set  2.  The 
resulting  values  can  be  seen  in  Table  V.  All  four  sets  of  code  converge  to  steady  state 
solutions,  with  larger  differences  than  what  were  observed  in  case  2,  but  had  only 
minor  differences  in  n'max,  vr 'min,  umax,  and  umin.  The  most  predominate  differences 
were  observed  in  the  wmax  and  wmin  between  x-z  coordinates  and  x-oy  coordinates, 
where  x-oy  had  greater  magnitude  values  for  both  state  variables.  There  was  also  a 
noticeable  difference  in  6'max  and  0'min  between  x-z  coordinates  and  x-oy  coordinates 
for  both  sets.  Consistent  with  the  previous  two  cases,  the  CPU  time  was  slower  for 
both  x-oy  model  simulations  (  52.06  %  increase  in  time  for  set  1  and  15.12  %  increase 
in  time  for  set  2)  relative  to  the  x-z  simulations.  This  difference  is  again  being  caused 
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by  the  additional  x-az  data  being  passed  between  program  functions  and  additional 
calculations  being  performed  (  0(5  *  [NT]  *  [NE]  *  [Mat]2  *  [8  *  (N  +  1)  +  26])).  For 
this  case,  there  were  1.8700xl016  additional  operations. 


Figure  12.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Resulting  horizontal  velocity 
(a)  and  vertical  velocity  (b)  using  both  modified,  (ii),  and  unmodified,  (i),  CGI  source 
codes  [1]  after  5  h  for  the  resolution  of  360  m  (in  x)  and  300  m  (in  z)  (colored  lines) 
plotted  with  the  analytic  solution  (black  lines).  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.025  to  0.025  ms-1  with  an  interval  of  0.005  ms-1, 
(a),  and  -0.005  to  0.005  ms-1  with  an  interval  of  0.0005  ms-1,  (b). 


60 


ii) 


Figure  13.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Resulting  horizontal  velocity 
(a)  and  vertical  velocity  (b)  using  both  modified,  (ii),  and  unmodified,  (i),  CG2  source 
codes  [1]  after  5  h  for  the  resolution  of  360  m  (in  x)  and  300  m  (in  z)  (colored  lines) 
plotted  with  the  analytic  solution  (black  lines).  All  cases  were  run  using  10th  order 
polynomials,  with  contours  from  -0.025  to  0.025  ms”1  with  an  interval  of  0.005  ms”1, 
(a),  and  -0.005  to  0.005  ms”1  with  an  interval  of  0.0005  ms”1,  (b). 


To  further  evaluate  the  performance  of  the  four  model  runs,  the  Root  Mean 
Squared  Error  (RMSE)  was  determined  for  each  run  in  relation  to  the  analytic  solution 
(see  Table  VI)  for  the  four  variables:  n,  u,  w ,  and  9.  Once  again,  with  the  domain 
of  interest  defined  and  the  numerical  and  analytic  solutions  on  matching  grids  a 
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Table  V.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Comparison  of  modified  and 
unmodified  CGI  and  CG2  using  360  m  (in  x)  and  300  m  (in  z)  resolution  after  5  h. 


Model 

CGI  x-z 

CGI  x-crz 

CG2  x-z 

CG2  x-az 

7 t' 

11 max 

0.2281x  10-6 

0.2059xl0-s 

0.2276xl0-s 

0.2045x  10-6 

7T/  • 
mm 

-0.2631X10-6 

-0.2628X10-6 

-0.2599x  10-6 

-0.2654x  10-6 

umax  (ms-1) 

0.8257x  10-2 

0.8194xl0-2 

0.8229X10-2 

0.8191x  10-2 

Umin  (ms-1) 

-0.7163X10-2 

-0.6910xl0-2 

-0.7120x  10-2 

-0.6903x  10-2 

Wmax  (ms-1) 

0.6035X10-2 

0.6696xl0-2 

0.6035X10-2 

0.6694x  10-2 

Wmin  (ms-1) 

-0.6037xl0-2 

-0.7989xl0-2 

-0.6037xl0-2 

-0.7988x  10-2 

0'max  (K) 

0.2962x  10-2 

0.4275X10-2 

0.2969X10-2 

0.4262x  10-2 

0'min  (K) 

-0.2846X10-2 

-0.3987xl0-2 

-0.2846x  10-2 

-0.3991x  10-2 

CPU  time  (s) 

0.6166X104 

0.9376x  104 

0.1720x10s 

0. 1980x10s 

bootstrap  (random  sampling)  method,  using  10,000  iterations,  was  used  to  calculate 
the  95  %  Cl  [11].  The  resulting  confidence  intervals  again  indicate  that  the  differences 
between  the  state  variable  RMS  errors  are  significant,  since  the  ranges  do  not  overlap. 
Both  CGI  x-z  and  CG2  x-z  exhibited  a  smaller  error  for  all  four  variables  of  interest, 
with  the  most  significant  difference  in  RMSE  values  for  w  and  6.  Similar  to  the 
previous  case,  taking  a  closer  look  at  the  w  RMSE  with  respect  to  wmax  shows  that 
CGI  x-z  has  a  2.17%  error  as  compared  to  20.25%  for  CGI  x-oz  and  that  CG2  x-z 
has  a  2.17%  error  as  compared  to  20.26%  for  CG2  x-az,  which  is  even  larger  than 
observed  in  the  linear  hydrostatic  case. 

Similar  to  the  previous  case,  in  order  to  verify  that  a  steady-state  solution  for 
the  linear  non-hydrostatic  mountain  case  was  achieved,  the  momentum  flux  was  again 
derived  using  Eq.  (5.1).  [12]  For  this  case,  the  analytic  non-hydrostatic  momentum 
flux  is  given  by: 


mNH(z )  =  —  0A57mH  (z) 


where  mNH  denotes  the  analytic  non- hydrostatic  momentum  flux  and  rnH  denotes 
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the  analytic  hydrostatic  momentum  flux  (Eq.  (5.2)).  [13]  The  momentum  flux  was 
then  normalized  by  m(z)/mNH(z),  and  will  be  the  value  discussed  in  this  thesis.  In 
Figure  14,  the  normalized  momentum  flux  for  all  four  model  runs,  using  360  m  (in  x) 
and  300  m  (in  z)  resolution,  was  plotted  for  1  h,  2  h,  3  h,  4  h,  and  5  h.  The  resulting 
momentum  flux  for  this  case  resembled  the  same  pattern  as  seen  in  case  2.  Figure  14 
indicates  that  all  four  models  converge  to  a  steady  state  solution  after  5  h.  Of  the 
four  model  simulations,  the  CGI  x-z  and  CG2  x-z  models  yielded  results  that  were 
far  better  than  the  x-az  models,  with  values  between  0.95  and  1.0  versus  between 
0.85  and  1.18  as  seen  for  the  x-crz  runs. 

2.  Conclusions 

All  four  model  runs  for  the  linear  non-hydrostatic  mountain  case  did  develop 
a  steady-state  mountain  wave  over  a  single  peak  indicating  that  the  model  dynamics 
are  resolving  the  scenario.  Additionally,  the  models  converged  to  nearly  identical 
patterns,  but  with  varying  oscillation  intensities.  Similar  to  the  results  seen  in  case  2, 
for  both  set  1  and  set  2  of  the  Navier-Stokes  equations,  the  x-az  coordinates  performed 
worse  than  their  x-z  counterparts  (reflected  in  the  RMSE,  the  normalized  momentum 
flux  ,  and  most  pronounced  in  w  with  an  approximately  nine  times  larger  RMSE 
with  respect  to  wmax  )  and  did  show  a  noticeable  increase  in  computational  expense. 
These  results  further  indicate  that  when  the  model  dynamics  are  resolved  using  purely 
explicit  time  integration  methods  that  there  is  no  significant  benefit  in  using  x-cr2  over 
x-z  and  that  there  is  a  significant  degradation,  but  the  results  are  promising  since  the 
models  are  converging  to  a  fairly  representative  steady-state  solution.  Although  the 
filters  and  boundary  conditions  are  designed  to  be  independent  of  z  and  az,  it  cannot 
be  ruled  out  that  a  special  treatment  could  be  explored.  With  more  research  into  the 
filters  and  boundary  conditions,  the  degradation  could  be  made  minimal  enough  that 
x-crz  is  worth  using  with  semi-implicit  methods. 
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Table  VI.  Case  3:  Linear  Non-Hydrostatic  Mountain  RMSE.  Root-mean-squared  er¬ 
rors  for  the  four  primary  state  variables,  for  both  the  modified  and  unmodified  codes, 
after  5  h  using  360  m  (  in  x)  and  300  m  (in  z)  resolution  and  10th  order  polynomials. 


Model 

7T 

95%  Conhdence  Interval 

CGI  x-z 

1.5939x  10-8 

+  1.0987xl0-lu 

-  1.1296xlO-10 

CGI  x-az 

2.2354x  10-8 

+  1.8145xlO-10 

-  1.8058xl0-1° 

CG2  x-z 

1.6391x  10-8 

+  1.1292xlO-10 

-  1.1720xl0-1° 

CG2  x-az 

2.2931x  10-8 

+  1.8913xlO-10 

-  1.8253xlO-10 

Model 

u  (ms-1) 

95%  Conhdence  Interval 

CGI  x-z 

4.6827x  10-4 

+  3.5124xl0-ti 

-  3.5127xlO-B 

CGI  x-az 

5.3962x  10-4 

+  4.3611X10-6 

-  4.3393x  10-6 

CG2  x-z 

4.8941  xlO-4 

+  3.6237xl0-6 

-  3.6915x  10-6 

CG2  x-az 

5.3838x  10-4 

+  4.3711X10-6 

-  4.3818x  10-6 

Model 

w  (ms-1) 

95%  Conhdence  Interval 

CGI  x-z 

1.3125x  10-4 

+  1.2912X10-6 

-  1.3101x  10-6 

CGI  x-az 

1.6180x  10-3 

+  1.2733X10-5 

-  1.2822x  10-5 

CG2  x-z 

1.3105x  10-4 

+  1.3205X10-6 

-  1.3154x  10-6 

CG2  x-az 

1.6180x  10-3 

+  1.2727xl0-5 

-  1.2691x  10-5 

Model 

d(K) 

95%  Conhdence  Interval 

CGI  x-z 

1.6683x  10-4 

+  1.3953X10-6 

-  1.4053x  10-6 

CGI  x-az 

5.3367xl0-4 

+  4.7779X10-6 

-  4.7821x  10-6 

CG2  x-z 

1.7285x  10-4 

+  1.4982X10-6 

-  1.4429x  10-6 

CG2  x-az 

5.3364x  10-4 

+  4.7514xl0-6 

-  4.7686x  10-6 
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Figure  14.  Case  3:  Linear  Non-Hydrostatic  Mountain.  Normalized  momentum  flux 
for  the  resolution  of  360  m  (in  x)  and  300  m  (in  z),  at  times  1  h,  2  h,  3  h,  4  h,  and 
5  h  for  the  four  model  runs:  (i)  CGI  x-z,  (ii)  CGI  x-az,  (iii)  CG2  x-z,  and  (iv)  CG2 
x-oy. 
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VI.  CONCLUSIONS  AND 
RECOMMENDATIONS 

The  results  from  the  three  test  cases  yielded  promising  results.  The  first  of 
which  was  that  the  x-a~  coordinates  functioned  properly  and  did  in  fact  reduce  to  x-z 
coordinates  in  the  absence  of  terrain,  as  demonstrated  in  case  1.  With  the  introduc¬ 
tion  of  terrain,  all  four  models  converged  to  nearly  identical  steady  state  patterns  for 
both  case  2  and  case  3  respectively,  but  had  varying  oscillation  intensities  between 
coordinate  systems.  Additionally,  both  case  2  and  case  3  x-crz  model  runs  did  de¬ 
velop  a  steady-state  mountain  wave  over  a  single  peak  (reflected  in  the  normalized 
momentum  flux)  indicating  that  the  model  dynamics  are  resolving  the  scenarios  for 
a  linear  hydrostatic  mountain  and  a  linear  non-hydrostatic  mountain.  In  general,  the 
CGI  x-crz  and  CG2  x-<r2  models  performed  worse  than  their  x-z  counterparts.  The 
x-crz  models  had  higher  RMS  errors  (as  large  as  nine  times  greater  RMSE  values  with 
respect  to  associated  maximum  vertical  velocities),  which  were  observed  predomi¬ 
nantly  in  intensity  values  and  not  in  placement  of  steady  state  features.  All  three 
cases  for  x-az  also  showed  a  noticeable  increase  in  computational  expense,  due  to 
the  additional  calculations  and  variables  required  by  the  coordinate  transformation. 
These  results  indicate  that  though  x-crz  coordinates  are  not  as  accurate  or  efficient 
as  x-z  and  that  there  is  a  significant  degradation,  with  further  fine-tuning  of  the 
model  environment  these  issues  could  be  made  minimal  enough  that  it  justifies  their 
use  in  semi-implicit  methods,  especially  in  the  vertical,  as  is  already  done  by  various 
operational  mesoscales  models. 

Past  research  has  demonstrated  the  utility  of  using  Continuous  Galerkin  meth¬ 
ods,  in  a  x-z  framework,  for  resolving  computational  fluid  dynamics  using  both  fully 
explicit  and  semi-implicit  time  integration  methods.  This  thesis  brought  Giraldo’s 
2-D  (x-z  slice)  mesoscale  Non-Hydrostatic  model  one  step  closer  to  evaluating  and 
exploiting  the  full  strength  of  x-az  coordinates  by  transforming  the  Navier-Stokes 
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Equations  and  testing  their  impacts  using  fully  explicit  time  integration.  Now  that 
the  x-az  models  are  functional,  the  next  stage  of  development  will  be  to  implement  a 
semi-implicit  time  integration  method  in  the  vertical  (1-D)  and  compare  the  results 
against  the  2-D  semi-implicit  x-z  model.  Without  the  work  done  in  this  project  to 
transform  the  x-z  equations  to  x-az,  one  could  not  construct  semi- implicit  methods  in 
the  vertical  since  in  x-z,  the  terrain  is  coupled  to  box  coordinates  and  so  the  coordi¬ 
nates  cannot  be  decoupled.  Additionally,  the  mathematically  machinery  outlined  in 
this  thesis  can  be  used  to  transform  any  equation  set  to  any  other  coordinate  system. 

The  value  of  this  study  is  far  reaching  in  determining  the  usefulness  of  ap¬ 
plying  a  specific  coordinate  system  in  the  future  when  developing  meteorological  and 
oceanographic  models  for  the  U.S.  Naval  Research  Laboratory  (NRL)  by  constituents 
at  the  Naval  Postgraduate  School.  In  addition,  the  successful  conversion  of  the  non¬ 
hydrostatic  x-z  models  to  x-az  will  allow  for  the  straightforward  extension  of  these 
models  to  global  non-hydrostatic  models,  since  cr  will  represent  the  height  of  the 
model  and  x  will  then  represent  the  spherical  shell  at  each  value  of  cr. 
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APPENDIX.  COEFFICIENTS  FOR  RK35 

METHOD 


Coefficient  Value 


stage  1: 

®  0 

1 

A) 

0.377268915331368 

stage  2: 

ot  0 

0 

ft  i 

1 

Po 

0.377268915331368 

stage  3: 

«o 

0.355909775063327 

0 

02 

0.644090224936674 

Po 

0.242995220537396 

stage  4: 

«o 

0.367933791638137 

o  1 

0 

«2 

0 

ft3 

0.632066208361863 

Po 

0.238458932846290 

stage  5: 

fto 

0 

O'! 

0 

ft2 

0.237593836598569 

ft3 

0 

ft4 

0.762406163401431 

Po 

0.287632146308408 
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